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ABSTRACT 


This  report  presents  a  satellite  simulation  algorithm  written  for  the  IBM  PC  and  clones 
in  Microsoft  QuickBASIC  4.5.  The  program  can  either  generate  an  ephemeris  of  satellite 
ground-track  position  and  location  with  respect  to  a  fixed  ground  station,  or  determine  the 
position  and  time  the  satellite  is  above  the  horizon  of  a  specified  ground  station.  Provision 
is  made  for  the  program  to  handle  multiple  satellites  and  multiple  ground  stations.  Input 
is  either  by  way  of  the  computer  keyboard  or  by  separate  files  for  satellites  and  ground 
stations.  Output  is  written  both  on  the  screen  and  to  an  ASCII  output  file  which  can  be 
accessed  by  a  data  based  management  program. 
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I.  INTRODUCTION 


This  report  presents  a  satellite  simulation  algorithm  written  for  the  IBM  PC  and  clones 
in  Microsoft  QuickBASIC  4.5.  The  program  can  either  generate  an  ephemeris  of  satellite 
•  ground-track  position  and  location  with  respect  to  a  fixed  ground  station,  or  determine  the 

position  and  time  the  satellite  is  above  the  horizon  of  a  specified  ground  station.  Provision 
v  is  made  for  the  program  to  handle  multiple  satellites  and  multiple  ground  stations.  Input 

is  either  by  way  of  the  computer  keyboard  or  by  separate  files  for  satellites  and  ground 
stations.  Output  is  written  both  on  the  screen  and  to  an  ASCII  output  file  which  can  be 
accessed  by  a  data  based  management  program. 

Two  versions  of  the  satellite  kinematics  subroutine  pos  .update  are  presented.  These 
algorithms  were  provided  by  NAVSPASUR  [Ref.  1].  The  first  is  a  very  accurate  algorithm 
based  upon  a  paper  by  Brouwer  [Ref.  2]  and  is  contained  in  NAVSPASUR's  FORTRAN  pro¬ 
gram  called  SHOWALL.  A  listing  of  the  BASIC  version  of  this  algorithm  is  presented  in 
Appendix  E.  The  second  algorithm  is  a  much  simplified  approximation  to  the  Brouwer 
model.  Equations  for  the  simplified  model  are  contained  in  Appendix  A,  and  the  BASIC 
implemention  is  embedded  in  the  listing  of  the  full  simulation  program  in  Appendix  D. 

The  next  sections  contain  program  operating  instructions,  structure  of  the  input/ 
output  files,  and  suggestions  for  future  improvements.  The  computation  of  the  times  of 
culmination,  rise  and  set  are  in  Appendix  B,  and  the  computation  of  satellite  heading  is  in 
Appendix  C.  Readers  unfamiliar  with  the  nomenclature  of  satellite  orbital  elements  should 
consult  any  standard  work  such  as  Escobal  [Ref.  3,  Ch.  3]. 

II.  OPERATING  INSTRUCTIONS 

A.  Overview.  The  source  code  is  named  satsta.bas  and  the  executable  code  is  named 
satsta.exe.  To  run,  enter  SATSTA  at  the  DOS  prompt  and  press  the  *_j  key. 

The  first  menu  (Fig.  1)  requests  a  data  input  mode.  Pressing  1  or  2  will  require  the 
data  for  one  satellite  and  one  ground  station  to  be  input  from  the  keyboard.  Pressing  3  or  4 
assumes  that  two  data  files  are  accessible:  either  elements.dat  or  elements. nss  containing 
the  orbital  elements  for  one  or  more  satellites,  and  station.dat  containing  the  location  of 
one  or  more  ground  stations — the  structure  of  these  files  is  discussed  in  the  next  section. 
No  matter  which  of  these  four  options  is  chosen,  a  simulation  starting  time,  ending  time 
and  time-step  increment  will  be  required  input  from  the  keyboard  in  the  third  menu. 

The  second  menu  (Fig.  2)  requests  a  run  mode.  Option  1  creates  an  ephemeris  of 
satellite  geographic  position  and  satellite  position  with  respect  to  the  ground  station  for 
equally  spaced  times  for  the  entire  simulation  duration.  Option  2  computes  the  satellite 
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geographic  position  and  satellite  position  with  respect  to  the  ground  station  only  for  those 
times  during  the  simulation  duration  at  which  the  satellite  is  above  the  horizon  of  the 
ground  station. 

The  third  menu  (Fig.  3)  requests  the  simulation-times,  consisting  of  a  starting  time, 
an  ending  time,  and  a  time-step  increment  (At). 

The  fourth  and  fifth  menus  (Figs.  4  and  5,  respectively)  appear  only  if  the  first  menu 
(data  input  mode)  options  1  or  2  were  chosen.  Menu  four  requests  the  satellite  epoch  and 
elements,  and  menu  five  requests  the  station  coordinates. 

No  matter  which  options  are  chosen,  output  is  directed  to  the  screen  and  is  appended 
to  file  output.dat.  If  the  file  does  not  exist,  it  is  created.  This  ASCII  file  is  structured  so  that 
it  can  be  used  as  input  to  a  data  base  management  program,  if  desired.  The  structure  of 
this  file  is  discussed  in  a  later  section. 

B.  Input.  Menu  1  requests  the  selection  of  a  data  input  mode.  Depending  upon  the 
source  of  satellite  elements,  there  are  two  ways  in  which  these  data  can  be  input.  The 
first  is  using  standard  "epoch  of  date1'  elements  in  which  the  time  of  perigee  passage  is  the 
epoch.  The  second  is  the  NAVSPASUR  "One-Line  Charlie”  elements  in  which  the  longitude 
(or  right  ascension)  of  the  ascending  node  has  been  referenced  to  0000  hours,  1  January 
19S5.  The  "One-Line  Charlie”  option  has  been  added  for  those  users  who  may  obtain 
satellite  elements  from  NAVSPASUR.  An  example  of  “epoch  of  date”  input  is  shown  in 
Fig.  3. 


Select  data  input  mode  1,  2,  3  or  4: 

Keyboard  input  lor  satellite  ft  ground  station: 

1.  Epoch  of  date. 

2.  NAVSPASUR  Reference  Date. 

Data  file  input  for  satellite  A  ground  station: 

3.  Epoch  of  date. 

4.  NAVSPASUR  One-Line  Charlie. 

Figure  1.  Menu  1 — Select  data  input  mode. 


If  data  file  input  Option  3  is  chosen  from  Menu  1,  two  ASCII  data  files,  elements.dat  and 
station .dat,  must  have  been  previously  created.  The  elements.dat  file  must  contain  a  single 
line  of  for  each  satellite  of  interest.  Each  line  must  contain  10  entries,  each  enclosed  in 
double  quotations  and  separated  by  a  blank  space  (Fig.  6).  The  entries  are:  (1)  a  satellite 
identification  consisting  of  at  most  five  alphanumeric  characters;  (2)  the  epoch  date  (date 
of  perigee  passage)  in  the  form  yyyy.mmdd,  where  yyyy  is  the  year,  mm  is  the  month,  and 
dd  is  the  day;  (3)  the  epoch  time  (time  of  perigee  passage)  in  the  form  hh.mmss,  where  hh 
is  the  hour,  mm  is  the  minutes,  and  ss  is  the  seconds;  (4)  the  orbital  eccentricity;  (5)  the 
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longitude  (or  right  ascension)  of  the  ascending  node  in  degrees;  (6)  the  orbital  inclination 
in  degrees;  (7)  the  argument  of  perigee  in  degrees;  (8)  the  mean  anomaly  in  degrees; 
(9)  the  mean  motion  in  revolutions  per  day;  and  (10)  the  orbital  decay  in  revolutions  per 
day-squared. 

If  data  file  input  Option  4  is  chosen  from  Menu  1,  two  ASCII  data  files,  elements. nss 
and  station.dat,  must  have  been  previously  created.  The  elements. mss  file  must  contain  a 
single  line  of  for  each  satellite  of  interest  in  the  NAVSPASUR  format.  Each  line  must  contain 
10  entries  starting  in  column  1  and  ending  in  column  65  (Fig.  7).  The  entries  are:  (1)  a 
satellite  identification  consisting  of  at  most  five  alphanumeric  characters  (columns  1-5); 
(2)  the  mean  anomaly  in  revolutions  (columns  6-13);  (3)  the  mean  motion  in  radians  per 
Herg1  (columns  14-21);  (4)  the  orbital  decay  in  radians  per  Herg-squared  (columns  22- 
27);  (5)  the  orbital  eccentricity  (columns  28-35);  (6)  the  argument  of  perigee  in  revolutions 
(columns  36-43);  (7)  the  longitude  (or  right  ascension)  of  the  ascending  node  in  revolutions 
(columns  44-51):  (S)  the  orbital  inclination  in  revolutions  (columns  52-59);  and  (9)  the 
epoch  date  (date  of  perigee  passage)  in  the  form  yymmdd,  where  yyyy  is  the  year,  mm  is  the 
month,  and  dd  is  the  day  (columns  60-65).  A  decimal  point  must  be  present  in  columns  6, 
14.  22.  2S.  3G.  44  and  52. 

Select  run  mode  1  or  2 : 

1.  Create  an  epheoeris. 

2.  Find  timea  satellite  is  above  the  horizon. 

Figure  2.  Menu  2 — Select  run  mode. 

The  station.dat  file  must  contain  a  single  line  of  for  each  ground  station  or  geographical 
location  of  interest.  Each  line  must  contain  5  entries,  each  enclosed  in  double  quotations 
and  separated  by  a  blank  space  (Fig.  8).  The  entries  are;  (1)  a  station  identification 
consisting  of  at  most  three  alphanumeric  characters;  (2)  the  station  latitude  in  the  form 
dd.mmss  (minus  if  south),  where  dd  is  degrees,  mm  is  minutes,  and  ss  is  seconds;  (3)  the  sta¬ 
tion  longitude  in  the  form  ddd  .mmss  (minus  if  west);  (4)  the  station  altitude  in  feet  (minus 
if  below  sea  level);  and  (5)  the  name  of  the  station  or  station  location  (alphanumeric). 

If  a  data  file  input  option  is  chosen,  then  SATSTA  will  open  the  two  files  described 
above  and  compute  either  an  ephemeris  or  the  times  that  each  satellite  is  above  the  horizon 
(Option  1  or  2  of  the  second  menu)  for  each  ground  station.  Otherwise,  keyboard  input  for 
one  satellite  and  one  ground  station  will  be  requested  from  the  fourth  and  fifth  menues, 
respectively. 

'The  Herg — a  unit  of  time  named  after  the  astronomer  Paul  Herget — is  defined  as  the  period,  divided  by 
2tr,  of  a  fictitious  satellite  at  zero  altitude  over  a  smooth  spherical  earth  with  no  atmosphere. 
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Simulation-Tin*  Data: 

Starting  dat*  (yyyy.madd)  ?  1983.0201 

Starting  tim*  (hh.mm**)  ?  00.00 

Ending  dat*  (yyyy.madd)  7  1983.0202 

Ending  tim*  (hh.ma**)  7  00.00 

Tim*  increment  (minut**)  7  01 

Press  any  k«y  to  continu*. 

Figure  3.  Menu  3 — Select  simulation  times. 

Five  pieces  of  information  regarding  the  simulation  are  requested  from  the  third  menu 
(Fig.  3):  (1)  the  simulation  starting  date  in  the  form  yyyy  .mmdd,  (2)  the  simulation  starting 
time  in  the  form  hh.mmss;  (3)  the  simulation  ending  date  in  the  form  yyyy.mmdd;  (4)  the 
simulation  ending  time  in  the  form  hh.mmss;  and  (5)  the  simulation  time-step  increment 
in  minutes.  The  screen  will  display  sample  data  for  each  entry  in  the  appropriate  format. 
If  any  sample  entry  is  to  be  changed,  it  should  be  completely  overwritten  in  order  for  it 
to  be  read  in  properly.  When  each  entry  is  correct,  enter  it  into  the  computer  by  pressing 
the  <_i  key. 

If  the  keyboard  entry  option  was  chosen  from  the  first  menu,  then  two  more  menues 
will  appear,  otherwise  the  simulation  will  commence  with  input  being  obtained  from  the 
two  data  files. 

Satellite  Data: 

Epoch  date  (yyyy.mmdd)  ?  1983.0201 

Epoch  time  (hh.mmss)  7  00.00 

Input  mean  «I«m*nt*  of  epoch: 

Satellite  ID  Number  ?  11111 

Eccentricity  7  0.0005545 

Long,  ascending  nod*  (d*g)  7  272.43497 

Inclination  (dag)  ?  65.06057 

Arg.  of  perigee  (deg)  ?  295.41470 

Mean  anomaly  (dag)  7  258.10682 

Mean  motion  (r*v/day)  ?  15.44194 

Decay  (r*v/day*2)  ?  0 

Press  any  key  to  continue. 

Figure  4.  Menu  4 — Input  satellite  data. 

Menu  4  (Fig.  4)  requires  ten  separate  entries:  (1)  the  epoch  date  (date  of  perigee 
passage)  in  the  form  yyyy.mmdd,  where  yyyy  is  the  year,  mm  is  the  month,  and  dd  is  the 
day;  (2)  the  epoch  time  (time  of  perigee  passage)  in  the  form  hh.mmss,  where  hh  is  the 
hour,  mm  is  the  minutes,  and  ss  is  the  seconds;  (3)  a  satellite  identification  consisting  of  at 
most  five  alphanumeric  characters;  (4)  the  orbital  eccentricity;  (5)  the  longitude  (or  right 
ascension)  of  the  ascending  node  in  degrees;  (6)  the  orbited  inclination  in  degrees:  (7)  the 
argument  of  perigee  in  degrees;  (8)  the  mean  anomaly  in  degrees;  (9)  the  mean  motion  in 
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revolutions  per  day;  and  (10)  the  orbital  decay  in  revolutions  per  day-squared.  The  screen 
will  display  sample  data  for  each  entry  in  the  appropriate  format.  Again,  if  any  sample 
entry  is  to  be  changed,  it  should  be  completely  overwritten  in  order  for  it  to  be  read  in 
properly.  When  each  entry  is  correct,  enter  it  into  the  computer  by  pressing  the  «_j  key. 

Station  Coordinates: 

Station  Identification  Number  ?  001 

Station  Location  or  Name  ?  Daisy,  Tenn. 

Station  latitude  (dd.mmss.  South  minus)  ?  35.12 

Station  longitude  (ddd.mmss.  Vest  minus)  ?  -85.12 

Station  altitude  (feet)  ?  500.0 

Press  any  hey  to  continue. 

Figure  5.  Menu  5 — Input  ground  station  data. 

Menu  5  (Fig.  5)  requires  five  separate  entries:  (1)  a  station  identification  consisting  of 
at  most  three  alphanumei'c  characters;  (2)  a  station  location  or  name:  (3)  the  station  lati¬ 
tude  in  the  form  dd.mmss  (minus  if  south);  (4)  the  station  longitude  in  the  form  ddd.mmss 
(minus  if  west):  and  (5)  the  station  altitude  in  feet  (minus  if  below  sea  level). 

C.  Output.  Output  is  written  to  the  screen  and  is  appended  to  the  data  file  output.dat 
if  it  exists,  otherwise  file  output.dat  is  created.  A  sample  from  the  output.dat  file  is  shown 
in  Fig.  9 — the  structure  of  this  file  is  discussed  later. 

Sample  output  from  the  “create  an  ephemeris”  option  screen  is  shown  in  Fig.  10. 
The  apogee/perigee  output  is  approximately  valid  only  for  the  first  orbit.  In  particular, 
the  longitude  is  not  corrected  for  the  earth’s  rotation  during  the  current  orbit.  Also,  the 
ascending  node  and  argument  of  perigee  axe  not  corrected  for  their  rates  of  change  during 
the  current  orbit.  It  should  be  further  noted  that  for  nearly  circular  orbits,  the  height  of 
the  satellite  will  probably  get  lower  than  the  perigee  value  and  higher  than  the  apogee 
value  at  locations  other  than  the  perigee  and  apogee  positions — this  is  an  effect  of  the 
earth's  oblateness. 

The  main  body  of  the  ephemeris  output  consists  of  date,  time,  satellite  latitude, 
longitude  and  height,  followed  by  the  elevation  angle,  azimuth,  and  range  (in  kilometers) 
of  the  satellite  from  the  station.  The  satellite  “look  angle”  is  the  angle  between  the  satellite 
“ground-zero”  point  (usually  called  the  satellite  “sub-point”)  and  the  station.  The  last 
column  of  output  is  the  satellite  heading. 

A  sample  of  the  output  from  the  “times  satellite  is  above  the  horizon”  option  screen 
for  a  time-step  (At)  of  1  minute  is  shown  in  Fig.  11.  First  the  time  of  culmination  for 
the  current  orbit  is  computed  and  saved.  The  culmination  time  is  then  incremented  and 
decremented  by  At  until  the  satellite  is  below  the  horizon.  The  time  of  rise  and  set  is 
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Figure  9.  Sample  output.dat  file  (some  output  expunged).  Epoch  of  date. 


Satellite:  11111 


Station:  001  Daisy,  Tann. 


Apogee: 

Haight 

Latitude 

Longituda 

Perigee : 
Height 
Latitude 
Longitude 


442.8  km. 

-55.1540 

100.2099 


450.4  km. 

55.1539 

280.2099 


Tua  1  Fab  1983 


year 

mo 

da 

hr 

win 

sec  lat 

long 

height 
(km)  elev 

azim 

range 

look 

angle 

head 

1983 

2 

1 

0 

0 

0.0  -12.26 

327.56 

434.0  -31.60 

123.30 

7437 

52.90 

157.82 

1983 

2 

1 

0 

10 

0.0  -45.84 

347.68 

441.5  -50.59 

136.39 

10382 

36.56 

145.02 

1983 

2 

1 

0 

20 

0.0  -65.19 

48.39 

446.4  -69.07 

151.85 

12354 

19.68 

88.46 

1983 

2 

1 

0 

30 

0.0  -44.71 

106.99 

439.1  -83.25 

221.39 

13089 

6.45 

34.07 

1983 

2 

1 

0 

40 

0.0  -10.95 

126.55 

430.7  -70.51 

302.57 

12481 

18.11 

22.04 

1983 

2 

1 

0 

50 

0.0  24.08 

141.07 

435.9  -51.76 

316.67 

10574 

35.26 

24.38 

1983 

2 

1 

1 

0 

0.0  55.55 

168.97 

449.2  -32.06 

324.49 

7561 

52.31 

46.47 

Figure  10.  Output  from  ephemeris  option. 

Epoch  of  date. 

Satellite : 

11111 

Station 

:  001 

Daisy,  Tenn. 

year 

mo 

da 

hr 

mn 

sec 

lat 

long 

height 

(km) 

elev 

azim 

range 

look 

angle 

head 

1983 

2 

1 

1 

13 

19.4 

54.37 

263.10 

450.1 

0.00 

340.48 

2436 

69.23 

135.43 

1983 

2 

1 

1 

13 

42.1 

53.35 

264.76 

449.7 

1.39 

341.68 

2285 

69.19 

136.92 

1983 

2 

1 

1 

14 

42.1 

50.53 

268.75 

448.7 

5.55 

345.74 

1892 

68.57 

140.46 

1983 

2 

1 

1 

15 

42.1 

47.58 

272.27 

447.6 

10.78 

351.87 

1514 

66.79 

143.50 

1983 

2 

1 

1 

16 

42.1 

44.51 

275.37 

446.4 

17.70 

2.15 

1167 

63.08 

146.09 

1983 

2 

1 

1 

17 

42.1 

41.35 

278.14 

445.1 

26.48 

21.39 

892 

56.95 

148.32 

1983 

2 

1 

1 

18 

42.1 

38.12 

280.62 

443.8 

32.24 

55.33 

773 

52.35 

150.23 

1983 

2 

1 

1 

19 

42.1 

34.83 

282.88 

442.6 

27.03 

90.21 

875 

56.40 

151.86 

1983 

2 

1 

1 

20 

42.1 

31.49 

284 . 94 

441.3 

18.10 

110.45 

1141 

62.69 

153.26 

1983 

2 

1 

1 

21 

42.1 

28.11 

286.84 

440.1 

10.97 

121.17 

1484 

66.60 

154.46 

1983 

2 

1 

1 

22 

42.1 

24.69 

288.61 

438.9 

5.61 

127.48 

1862 

68.51 

155.46 

1983 

2 

1 

1 

23 

42.1 

21.25 

290.27 

437.9 

1.36 

131.58 

2255 

69.20 

156.31 

1983 

2 

1 

1 

24 

3.8 

19.99 

290.85 

437.5 

-0.00 

132.74 

2400 

69.25 

156.58 

Figure  11.  Output  from  “times  above  horizon”  option.  A t  —  1  minute.  Epoch  of  date. 


iterated  until  the  time  at  which  the  satellite  is  on  the  horizon.  All  of  these  times  and 
positions  are  then  output  to  the  screen.  Note  that  the  information  printed  out  in  this 
option,  with  the  exception  of  the  apogee/perigee  positions,  is  identical  to  that  printed 
out  in  the  “ephemeris”  option.  By  making  A t  large,  say  10  minutes,  in  the  “times  above 
horizon”  option,  then  only  the  times  of  satellite  rise,  culmination,  and  set  axe  computed 
and  output  (see  Fig.  12). 

Each  line  in  the  output.dat  file  consists  of  the  station  identifier  and  the  satellite  iden¬ 
tifier  followed  by  the  same  data,  and  in  the  same  order,  as  that  printed  to  the  screen.  This 
can  be  seen  by  comparing  the  first  six  lines  of  Fig.  12  with  the  same  lines  in  Fig.  9. 
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Satellite:  11111  Station:  001  Daisy,  Tens. 


year 

mo 

da 

hr 

mn 

sec 

lat 

long 

height 

(km) 

elev 

azim 

look 

range  angle 

head 

1983 

2 

1 

1 

13 

19.4 

54.37 

263.10 

450.1 

0.00 

340.48 

2436 

69.23 

135.43 

1983 

2 

1 

1 

18 

42.1 

38.12 

280.62 

443.8 

32.24 

SS.33 

773 

62.35 

150.23 

1983 

2 

1 

1 

24 

3.8 

19.99 

290.85 

437.  S 

-0.00 

132.74 

2400 

69.25 

156.58 

1983 

2 

1 

2 

49 

40.9 

45.38 

250.95 

446.7 

0.00 

306.66 

2428 

69.23 

145.41 

1983 

2 

1 

2 

54 

12.9 

30.63 

261.84 

441.0 

11.96 

261.45 

1429 

66.15 

153.58 

1983 

2 

1 

2 

58 

44.7 

15.05 

269.44 

436.2 

-0.00 

195.12 

2394 

69.26 

157.45 

1983 

2 

1 

14 

47 

13.6 

16.76 

285.33 

433.6 

0.00 

150.29 

2387 

69.32 

22.83 

1983 

2 

1 

14 

50 

58.0 

29.64 

291.62 

437.9 

6.44 

106.44 

1792 

68.35 

26.06 

1983 

2 

1 

14 

54 

43.3 

42.02 

300.04 

443.2 

-0.00 

63.10 

2419 

69.28 

32.12 

1983 

2 

1 

16 

20 

48.3 

17.88 

262.22 

433.9 

0.00 

216.51 

2389 

69.32 

23.02 

1983 

2 

1 

16 

26 

16.7 

36.50 

272.26 

440.8 

85.22 

303.76 

528 

32.33 

28.94 

1983 

2 

1 

16 

31 

46.6 

53.37 

289.30 

448.2 

-0.00 

24.72 

2431 

69.26 

43.11 

1983 

2 

1 

17 

59 

48.3 

37.38 

249.28 

441.2 

0.00 

283.77 

2414 

69.29 

29.38 

1983 

2 

1 

18 

2 

51.4 

47.05 

257.60 

445.5 

3.50 

317.85 

2066 

69.01 

36.02 

1983 

2 

1 

18 

5 

55.6 

55.73 

269.75 

449.2 

-0.00 

351.87 

2433 

69.26 

46.77 

Figure  12.  Output  from  "times  above  horizon”  option.  A t  =  10  minutes.  Epoch  of  date. 


Fig.  13  is  a  sample  output  obtained  using  the  NAVSPASUR  reference  time.  A  compar¬ 
ison  of  Fig.  10  and  Fig.  13  shows  that  the  satellite  latitude,  height,  and  heading  is  the 
same  but  the  longitude  differs  because  of  the  difference  of  the  earth’s  rotation  between  the 
two  different  reference  times.  The  satellite-station  relationships  differ  for  the  same  rea¬ 
son.  Figures  12  and  14  are  completely  dissimilar  because  of  the  differing  satellite-station 
relationships. 


Satellite: 

uni 

Station 

:  001 

Daisy,  Tenn. 

Apogee : 
Height  - 

442 

.8  km. 

Latitude  = 

-55 

.1540 

Longitude  = 

230 

.8483 

Perigee: 
Height  * 

450 

.4  km. 

Latitude  * 

65 

.1539 

Longitude  - 

60 

.8483 

Tue  1  Feb  1983 

height 

look 

year  mo  da  hr  mn 

sec  lat 

long 

(km)  elev 

azim 

range  angle 

head 

1983  2  1 

0  0 

0.0  -12.26 

98.20 

434.0  -77.96 

351.78 

12919 

11.09 

157.82 

1983  2  1 

0  10 

0.0  -45.84 

118.32 

441.5  -79.40 

232.53 

12965 

10.02 

145.02 

1983  2  1 

0  20 

0.0  -65.19 

179.03 

446.4  -60.98 

210.44 

11620 

27.15 

88.46 

1983  2  1 

0  30 

0.0  -44.71 

237.62 

439.1  -41.32 

205 . 69 

9023 

44.80 

34.07 

1983  2  1 

0  40 

0.0  -10.95 

257.18 

430.7  -20.38 

203.40 

5458 

61.42 

22.04 

1983  2  1 

0  50 

0.0  24.08 

271.70 

435.9  12.73 

194.92 

1373 

65.77 

24.38 

1983  2  1 

1  0 

0.0  55.55 

299.61 

449.2  -5.07 

32.07 

3062 

68.64 

46.47 

Figure  13.  Output  from  ephemeris  option,  navspasur  reference  time. 
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Satellite:  11111  Station:  001  Daisy,  Tenn. 


height  look 


year 

mo 

da 

hr 

mn 

sac 

lat 

long 

(km) 

alav 

azim 

ranga 

angla 

haad 

1983 

2 

1 

0 

47 

32.1 

15.54 

267.79 

433.2 

0.00 

199.85 

2385 

69.32 

22.63 

1983 

2 

1 

0 

52 

59.8 

34.24 

277.36 

439.8 

58.70 

112.56 

509 

29.01 

27.87 

1983 

2 

1 

0 

58 

28.6 

51.41 

293.03 

447.4 

-0.00 

33.16 

2429 

69.27 

40.56 

1983 

2 

1 

2 

24 

57.1 

29.86 

251.15 

438.0 

0.00 

262.00 

2405 

69.30 

26.14 

1983 

2 

1 

2 

29 

10.7 

43.72 

260.91 

443.9 

8.69 

312.74 

1641 

67.65 

33.31 

1983 

2 

1 

2 

33 

25.7 

55.99 

276.83 

449.3 

-0.00 

3.07 

2433 

69.26 

47.22 

1983 

2 

1 

8 

58 

59.6 

56.01 

272.80 

450.7 

0.00 

356.74 

2437 

69.23 

132.74 

1983 

2 

1 

9 

3 

16.0 

43.68 

288.80 

446.0 

8.83 

47.32 

1637 

67.55 

146.72 

1983 

2 

1 

9 

7 

31.7 

29.71 

298.62 

440.7 

-0.00 

98.40 

2412 

69.24 

153.92 

1983 

2 

1 

10 

33 

57.2 

51.41 

256.65 

449.0 

0.00 

326.74 

2433 

69.23 

139.44 

1983 

2 

1 

10 

39 

26.8 

34.21 

272.34 

442.3 

58.17 

247.39 

514 

29.48 

152.14 

1983 

2 

1 

10 

44 

55.7 

15.45 

281.93 

436.3 

-0.00 

160.33 

2394 

69.26 

157.39 

1983 

2 

1 

12 

13 

31.5 

31.47 

250.42 

441.3 

0.00 

266.88 

2414 

69.24 

153.27 

1983 

2 

1 

12 

14 

27.8 

28.30 

252.20 

440.2 

0.29 

256.89 

2378 

69.24 

154.39 

1983 

2 

1 

12 

15 

23.2 

25.15 

253.85 

439.1 

-0.00 

247.03 

2407 

69.25 

155.34 

Figure  14.  Output  from  "times  above  horizon”  option.  A t  =  10  minutes.  Navspasur  reference  time. 


III.  SOURCES 

As  mentioned  in  the  introduction,  the  high  resolution  and  simplified  satellite  kinematics 
models  were  obtained  from  NAVSPASUR  [Ref.  1].  The  high  resolution  model  is  based  upon 
a  paper  by  Brouwer  [Ref.  2j.  The  code  for  the  high  resolution  model  is  given  Appendix  E, 
while  the  code  for  the  simplified  model  is  contained  in  the  pos. update  subroutine  in 
Appendix  D.  Both  of  these  routines  require  a  solution  of  Kepler’s  equation.  A  non-iterative 
state-of-the-art  solution  of  Kepler’s  equation,  based  upon  a  paper  by  Mikkola  [Ref.  6]  is 
contained  in  subroutine  kepler  listed  in  Appendix  D. 

The  ground  track  determination  is  essentially  that  given  by  Escobal  [Ref.  3,  Ch.  10.3.2) 
but  with  a  minor  improvement  described  by  Shudde  [Ref.  7). 

The  method  of  computation  of  the  times  of  satellite  culmination,  rise,  and  set  are  given 
in  Appendix  B,  and  the  computation  of  the  satellite  heading  is  contained  in  Appendix  C. 
Additional  references  are  given  there. 

The  IAU  (1976)  System  of  Astronomical  Constants  [Ref.  10,  pg.  S7]  is  used  in  SATSTA 
as  well  as  the  J2000.0  reference  frame  [Ref.  10,  pg.  LI].  Satellite  positions  computed  by 
NAVSPASUR’s  SHOWaLL  [Ref.  1]  use  the  WGS  72  constants  and  the  J1900.0  reference  frame. 
For  this  reason,  the  agreement  between  SATSTA  and  SHOWALL  is  not  exact. 


IV.  DISCUSSION 

Although  the  position  updating  in  SATSTA  is  based  upon  the  satellite  kinematics  model 
of  Brouwer  [Ref.  2],  other  models  are  currently  in  use.  Most  notable  is  the  work  of  Hoots 
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and  Roehrich  [Ref.  5].  SATRAK  [Ref.  4]  is  a  FORTRAN  program  implementing  their  models 
and  is  available  from  the  Aerospace  Defense  Command,  Peterson  AFB.  Unfortunately,  the 
source  code  is  not  available  for  modification. 

The  program  presented  here  can  be  used  as  a  basis  for  a  more  sophisticated  model. 
The  input  scheme  could  be  improved  by  using  an  input  editor  to  catch  improper  types  of 
entry  and  to  allow  corrections  or  changes  to  be  made  to  previous  entries.  File  management 
programs  can  be  used  to  select  satellite  elements  and  station  parameters  from  master  files 
and  insert  them  into  “working”  files  for  program  input. 

Various  sensor  types  have  not  been  modelled  and  no  provision  has  been  made  for 
including  them  in  the  program.  Also,  satellite-satellite  visablity  has  not  been  modelled.  If 
such  features  are  desirable,  they  could  be  included  by  modifying  the  existing  program. 
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APPENDIX  A.  THE  SIMPLIFIED  KINEMATICS  ALGORITHM. 

The  computation  of  the  perturbed  satellite  orbit  about  the  earth  is  performed  using  a 
method  which  is  an  approximation  of  that  used  by  NAVSPASUR  [Ref.  1].  According  to 
NAVSPASUR,  this  approximation  will  yield  satellite  positons  with  a  random  error  of  about 
10-15  nautical  miles  for  a  20  day  interval  around  epoch. 

The  mean  elements  at  epoch,  which  must  be  input,  are: 

Eccentricity  eo 

Right  ascension  of  the  ascending  node  Q0 

Inclination  io 

Argument  of  perigee  wo 

Mean  anomaly  Mo 

Mean  motion  Mi,  and 

Decay  M2 

Using  the  mean  elements,  the  following  auxiliary  quantities  are  computed: 


aQ  =  Mj 


do  =  do 


1  + 


C(3cos2  t'o  —  1) 


12/3 


Cj  —  C 


Q  =  -2  C 


«2(i  -  eo)3/2 

5  cos2  t'o  —  1 

aoC1  -  eo)2  ’ 

COS  t’o 


d  =  — 


ao(l  —  eo)2 
4a0 


and 


3  Mi 


M2, 


( A — 1 ) 


(A-la) 


(A-2) 


(A— 3) 


(A-4) 


where  C  =  z/aJ2  is  related  to  the  oblateness  of  the  earth.  NAVSPASUR  states  the  following: 

For  optimal  accuracy,  ao  from  Equation  A-l  is  used  in  Equation  A-la  to  obtain  a  final 
value  for  ao.  This  latter  ao  is  then  used  throughout.  Where  program  size  is  a  major 
consideration  and  additional  errors  of  some  ten  miles  in  satellite  height  can  be  tolerated 
the  user  may  choose  to  leave  out  Equation  A-la  and/or  Equation  A-4. 

Using  the  mean  values  of  the  elements  at  epoch  To  and  the  results  of  the  above 
calculations,  the  mean  values  of  the  elements  for  any  time  t  are: 


M  =  Mo  +  M,(t  -  To)  +  M2(t  -  T0)2, 
e  =  eo, 

=  u>o  +  w(M  —  Mo), 
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Cl  =  Oo  +  Cl(M  —  Mo)  —  u>e(t  —  To), 
i  =  io,  and 
a  =  a0  +  a(t  -  T0), 

where  u>e  is  the  earth’s  sidereal  rotation  rate. 
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APPENDIX  B.  THE  TIME  OF  CULMINATION,  RISE  AND  SET  ALGO¬ 
RITHM. 

In  this  program  the  earth  is  considered  to  be  non-rotating,  i.e.,  the  coordinates  of  a  station 
are  given  by  latitude  and  longitude  with  the  longitude  unaffected  by  the  earth’s  rotation. 
Rotation  is  accomplished  by  subtracting  the  accrued  rotation  from  the  right  ascension  of 
the  ascending  node  of  the  satellite.  Although  this  departs  from  realism,  the  benefit  is  that 
the  rectangular  coordinates  of  the  station  need  to  be  calculated  only  once. 


Referring  to  Fig.  B-l,  culmination  of  the  satellite,  S,  with  respect  to  a  ground  station, 
P,  occurs  when  the  elevation  angle,  h,  of  the  satellite  is  a  maximum,  which  is  equivalent 
to  the  minimization  of  the  of  the  zenith  angle,  90°  —  h,  which  is  also  equivalent  to  the 
minimization  of  the  angle  rj  between  z,  the  unit  vector  normal  to  the  station  location,  and 
r,  the  vector  from  the  earth’s  center,  0,  to  the  satellite.  The  equation  for  rj  is 


z  •  r 


cos  n  =  - , 

r 


(B-l) 


where  z  is  given  by 


z  = 


cos  4>  cos  A 
cos  <i>  sin  A 
sin  <i> 


(B-2) 


and  4>  and  A  are  the  station  geodetic  latitude  and  longitude,  respectively.  The  satellite 
position  vector  r  is  obtained  from  [Ref.  3,  Equ.  (3.57)] 


r  =  xJP  +  yu/Q, 


(B— 3) 
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where  xw  and  yu  measure  the  satellite  position  in  the  orbital  plane.1  P  and  Q  axe  the 
first  two  columns  of  the  Euler  matrix  [Ref.  3,  Equs.  (3.43)  and  (3.44)] — they  transform  the 
in-plane  xw  and  yw  coordinates  to  the  equatorial  rectangular  system  and  are  given  by 


(B-4a) 


>x' 

cos  u>  cos  T  —  sin  u  sin  T  cos  i 

p  = 

Py 

= 

cos  w  sin  T  +  sin  u>  cos  T  cos  i 

Pt\ 

sin  u>  sin  i 

and 


' Qz 

—  sin  u)  cos  T  —  cos  u  sin  T  cos  t 

Q  — 

Qy 

= 

—  sin  u  sin  T  4-  cos  ui  cos  T  cos  t 

[q,\ 

cos  ui  sin  t 

(B-Ab) 


where  u>  is  the  argument  of  perigee  and  i  is  the  orbital  inclination;  F  =  —  ut{t  -  M, 

where  ft  is  the  longitude  of  the  ascending  node  at  the  time,  to,  of  the  last  element  update; 
u>e  is  the  earth’s  sidereal  rotation  rate  and  t  is  the  current  time. 

The  components  of  r  vary  with  time.  We  will  transform  the  equations  in  such  a 
way  that  will  replace  time  by  the  eccentric  anomaly  E  as  the  independent  variable.  The 
inclination  i  is  time  independent  and  the  elements  ft  and  u >  will  be  considered  to  be 
independent  of  time  over  the  period  of  a  single  orbit  even  though  they  axe  time  dependent 
through  perturbations.  The  components  of  r  which  will  be  considered  to  be  time  dependent 
are  then  xw,  yw  and  T. 

Write  Equation  (B-l)  as 

z  •  r  —  r  cos  tj  =  0. 

Take  the  derivative  with  respect  to  E  to  obtain, 

dr  dr  drt 

2'dE-dEC°S’>  +  rSm,’jE=0- 

To  find  the  value  of  E  which  minimizes  rj,  set  drj/dE  =  0.  Thus, 

dr  dr 

*-3£-sros,,=0’ 


or 


dr  l  dr 

z  •  — - — z  •  r  =  0. 

dE  r  dE 

Since  this  equation  cannot  be  solved  explicitly  for  E ,  set 


F(E)  =  z 


dr  1  dr 
dE  r  dE 


z  •  r 


(B-5) 


1  The  origin  of  the  system  is  at  the  focus  of  the  orbit.  xw  is  positive  in  the  direction  of  the  perifocus.  yw 
lies  in  the  orbital  plane  and  is  orthogonal  to  zw. 
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and  find  E  such  that  F(E )  =  0  by  the  Newton-Raphson  method.  The  Newton- Raphson 
procedure  also  requires  F'(E),  the  derivative  of  F(E)  with  respect  to  E: 


_  d2r  1  f  dr  \2  1  d2r  1  dr  dr 

^  ^  Z'dE2+r2\dEJ  r  dE2  Z  T  rdEZ'dE' 


The  first  and  second  derivatives  of  Equation  (B-3)  are 


—  -  dx“  -D  ,  dp  ,  dy*  n  ,  dCl 
dE  dEP  +  X“dE+  d£  Q  +  Vul  dE 


and 


cPr  _  L  0dxy,dP  ,  _  d2P  ,  d2yw^  ,  „dyw  dQ  ,  d2Q 

'  ~  ’  —  *  —  d-  ^ IjJ  i  no  d-  ,  ,-,0  Q  d-  2  .  )n  ■+  y^j  ■ 


dE2  dE 2  '  '  “  dE  dE 

where  [Ref.  3,  Equ.  (3.69)] 


dE2  dE2  ~  dE  dE 
xw  =  a(cos  E  —  e)  and 


’  dE 2  ’ 


yw  =  a\J  1  —  e2  sin  E, 


with  first  and  second  derivatives 


dxu  . 

—  =  -asiaE, 


dy , 


— ^  =  a\/l  —  e2  cosE, 
db 


and 


d2yu 


d£2 


=  —  ay/l  —  e 2  sin  E  =  — yw, 


d£2 


=  —a  cos  E. 


The  first  and  second  derivatives  of  Equations  (B-4)  are 


dPT 


dr 


=  ~p  — 

*  V  ml 


dE  *vdE' 


dPL  =  pdT 
dE  xdE ’ 


dPx 


dE 


=  0, 


dQ^  _  dr 

dE;  V,d£’ 


dQy  _ 


=  0  ^ 

dE  WzdE' 


dQz 


dE 


=  0, 


(B-6) 

(B~7a) 

(B-7b) 

(B-8a) 

(B-Sb) 

(B-Sc) 

(B-9a) 
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and 


d2P , 
d£2 
d2Qr 
dE2 


=  0, 


=  -<?* 


d2Qv  _ 


dE2 

d2Qz 

dE2 


=  -Q* 


=  o. 


Also  required  is  [Ref.  3,  Equ.  (3.67)] 


and 


r  =  a(l  —  e  cos  E) 


(B-9b) 


(B-10a) 


and  it’s  derivatives 


=  ae  sin  E ,  and 
dE 


d2r 

dE2 


=  ae  cos  E. 


Finally,  using  Kepler’s  equation  [Ref.  3,  Equ.  (3.79)], 


E  —  e  sin  E  _ 

t  =  - +  T, 

n 


(B-lOb) 


( B— 1 1) 


where  n  is  the  rate  of  orbital  motion  and  T  is  the  time  of  last  perifocal  passage,  the 
equation  for  T  can  be  written  as 


Then, 


and 


_  _  ( E  —  e  sin  .E  _  \ 

T  =  Q,  —  u>e  ( - — - +  T  —  to  1 

dr  ( 1  —  e  cos  E  \ 

—  =  -we  ^  ~  j  . 


(B-12a) 


(B-12b) 

(B-12c) 


d2r  _  — u;e(esin^) 
dE2  “  n  ' 

To  solve  Equation  (B-5),  an  initial  estimate  of  E  is  made  and  the  right-hand  side 
of  Equation  (B-13)  is  evaluated.  The  result,  the  left-hand  side  of  Equation  (B-13),  is  an 
improved  estimate  of  E, 

F(E) 


E  =  E- 


F'(E)' 


(B-13) 
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This  iteration  is  continued  until  \F{E)/F'(E)\  is  less  than  some  prescribed  amount.  The 
value  of  E  obtained  will  either  maximize  or  minimize  77.  The  desired  value  of  77  will  be  the 
one  for  which  cos  77  >  0. 

Once  the  eccentric  anomaly  for  culmination,  Ec ,  is  found,  the  time  of  culmination  is 
determined  directly  from  Equation  (B-ll).  Perturbing  Ec  by  plus  or  minus  about  10°,  it 
is  possible  to  obtain  good  starting  values  for  the  determination  of  the  satellite  set  and  rise 
times,  respectively. 

The  procedure  for  finding  the  rise  and  set  times  is  very  similar  to  finding  the  time 
of  culmination.  The  rise/set  equations  are  developed  by  Escobal  [Ref.  3,  Sec.  5.4],  but 
the  algorithm  presented  here  is  somewhat  simplified  since,  unlike  Escobal,  the  time  of 
culmination  is  available.  In  fact,  many  of  the  same  equations  can  be  used.  Again,  the 
eccentric  anomaly  will  be  used  as  the  independent  variable.  Referring  to  Fig.  B-l,  satellite 
rise  or  set  occurs  when  the  elevation  angle  h  is  equal  to  zero,  where  h  can  be  found  from 
the  equation 

sin  h  =  — -,  (B-14) 

P 

where 

p  =  r  -  R.  (B— 15) 

In  Equation  (B— 15),  R  is  the  station  coordinate  vector  and  p  is  the  slant  range  vector 
from  the  station  to  the  satellite.  For  the  spheroid  earth  [Ref.  3,  Equ.  (5.54)], 


R  = 


G\  cos  <j>  cos  A 
G\  cos  <j>  sin  A 
G2  sin  <j> 


(B-16) 


where  <f>  and  A  are  the  station  geodetic  latitude  and  longitude,  respectively.  G\  and  G2 
are  factors  which  account  for  the  earth’s  oblateness,  and  are  given  by 


G\  =  t  0<!  +  H  and 

yl  _  (2/  -  /2)sin2  <f> 

g2  =  —  - +  H , 

\/l  -(2/-/2)sin2<£ 


(B— 17) 


where 


at  =  earth’s  equatorial  radius, 

/  =  earth’s  flattening  factor,  and 
H  =  station  elevation  above  the  ellipsoid. 
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Write  Equation  (B-14)  as 


p  ■  z  =  psin  h. 


When  the  satellite  is  on  the  horizon,  h  =  0  and  sin  h  =  0.  Since  the  resulting  equation, 
p  ■  z  =  0,  cannot  be  solved  explicitly  for  E,  set 

G(E)  =  p-  z 


and  find  E  such  that  G(E)  =  0  by  the  Newton-Raphson  method.  Using  Equation  (B-15), 
write  G(E)  as 

G(E)  =  z  ■  (r  -  R).  (B-18) 

The  Newton-Raphson  procedure  also  requires  G'(E),  the  derivative  of  G(E)  with  respect 
to  E.  Since  the  earth  is  being  treated  as  non-rotating  with  the  sidereal  rotation  subtracted 
from  the  satellite’s  ascending  node,  R  is  independent  of  E,  hence 


G\E) 


(B— 19) 


In  Equations  (B-18)  and  (B-19),  z,  r,  dr/dE ,  and  the  additional  needed  derivatives  have 
all  been  given  previously  in  this  section. 

To  solve  Equation  (B-18),  and  initial  estimate  of  E  is  made  and  the  right-hand  side 
of  Equation  (B-20)  is  evaluated.  The  result,  the  left-hand  side  of  Equation  (B-20),  is  an 
improved  estimate  of  E , 


E  =  E  - 


G(E) 

G'{E)' 


(B-20) 


The  iteration  is  continued  until  \G{E) / G' (E)\  is  less  than  some  prescribed  amount.  The 
corresponding  time  of  rise  or  set  is  then  obtained  from  Equation  (B-ll). 


APPENDIX  C.  THE  SATELLITE  HEADING  COMPUTATION. 


A  satellite  can  be  located  by  a  position  vector  x  =  xji  +  x2 j  +  13k  in  a  rectangular 
coordinate  system  with  origin  at  the  earth’s  center.  In  matrix  form,  and  with  x3,  x2  and 
X3  expressed  in  spherical  coordinates, 


x  — 


r  cos  <f>  cos  A 
r  cos  <i>  sin  A 
r  sin  <p 


where  <f>  is  the  latitude,  A  is  the  longitude  at  the  point,  and  r  =  jx|  is  the  distance  of  the 
point  from  the  earth’s  center. 

Let  z  be  the  unit  vector  in  the  direction  of  r.  Then 


- 

Zi 

cos  <f>  cos  A 

22 

— 

cos  4>  sin  A 

23  _ 

si^  ‘ 

It  is  convenient  to  define  two  other  unit  vectors  in  the  plane  tangent  to  z.  These  are  local 
north  n  and  local  east  e,  defined  by 


where 


II 

o|”  0 

N  |_S_ 

II 

S 

—  sin  (j>  cos  A 

—  sin  <p  sin  A 

cos  <f> 

and 

6  |za| 

—  sin  A 
cos  A 

0 

dz 

Z,t,  =  d4> 

and 

dz 

Zx  =  ex' 

The  vectors  z,  n  and  e  provide  the  basis  for  a  right-handed  orthogonal  coordinate  system. 

If  a  satellite’s  velocity  vector  v  =  vii  +  u2j  +  V3IC  is  known  at  position  x  then  it  is 
easily  shown  [Ref.  8]  that 

n  •  v  =  cos  a  and 


e  •  v  =  sin  a, 

where  a  is  the  course  angle  or  heading.  FYom  these  relations,  a  is  found  to  be 


a  =  qatn(e  •  v,  n  •  v) 

where  qatn(y,x)  is  arctan(y/x)  corrected  to  the  proper  quadrant. 

Since  the  position  and  velocity  component’s  are  in  rectangular  coordinates,  it  will  be 
necessary  to  convert  them  to  spherical  coordinates  in  order  to  determine  the  components 
of  n  and  e.  However,  the  spherical  coordinate  step  can  be  circumvented  as  follows:  First, 


cos  a  =  n  •  v  =  —  uj  sin  <f>  cos  A  —  v2  sin  <j>  sin  A  +  t>3  cos  <j>. 
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Multiply  both  sides  by  cos  <f>  to  give 


cos  <f>  cos  a  =  —  v\  sin  <p  cos  <f>  cos  X  —  v2  sin  <f>  cos  <f>  sin  A  +  v3  cos2  <j> 

=  — sin^(^iyi  4-  z2v2)  +  V3  cos2  <t>  4-  Z3V3  sin^  —  231*3  sin<£ 
=  —  sin^(2iux  4-  z2v2  4-  23V3)  4-  v2  cos2  <f>  4-  23U3  sin<£ 

=  — 23(z  •  v)  4-  U3  cos2  <j>  +  v3  sin2  <j> 

=  v3  -  23(z  •  v). 


Similarly, 

so  that 


sin  a  =  e  •  v  =  —  i>i  sin  A  4-  v2  cos  A, 
cos  <j>  sin  a  =  v 2  cos  <f>  cos  A  —  vx  cos  <f>  sin  A 


Then, 


=  zxv2  -  z2v  1. 


a  =  qatn(cos  </>sin  a,  cos  ^  cos  a). 


Since  cos<t>  >  0  for  all  <f>  €  (— 7r/2,  tt/2),  the  cos<£  term  effectively  “cancels  out”  [Ref.  9] 
and  so  the  heading  is  determined  by 


a  =  qatn(sina:,  cos  a) 

=  qatn(2jV2  -  z2vx,v3  -  z3(z  ■  v)). 

This  equation  is  used  by  NAVSPASUR  [Ref.  1]  in  SHOWALL.  In  SATSTA,  the  equivalent  form, 

(x \v2-x2vx  x3(x-v)\ 

a  =  qatn  I  - - - ,  v3 - ~2 - I  , 

is  used. 
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APPENDIX  D.  THE  SATSTA  PROGRAM  LISTING. 


'SATSTA 

*  Satellite  Track.  03-03-88.  Rev.  09-26-89  A  0830. 

DECLARE  FUNCTION  acost  (r#) 

DECLARE  FUNCTION  asint  (xt) 

DECLARE  FUNCTION  atan2t  (yt,  x#) 

DECLARE  FUNCTION  decimalt  (vt) 

DECLARE  FUNCTION  dmodi  (xt,  mt) 

DECLARE  FUNCTION  dot#  (a#(),  bt()) 

DECLARE  FUNCTION  gmstt  (tut,  tmt) 

DECLARE  SUB  apo. peri. gee  (elemt(),  «l«aut()) 

DECLARE  SUB  culmnate  (zvt(),  sat ,  tOt,  a at ,  cosat&t,  alemutO) 

DECLARE  SUB  aular  (parit,  nodat,  inclt,  paat(),  valtQ) 

DECLARE  SUB  ground. track  (rt,  aat.dct,  aat.ltt,  sat.htt) 

DECLARE  SUB  hr.min.aac  (tima.of .dayt ,  hrX,  minX,  aact,  nX) 

DECLARE  SUB  jd. to. data  (jdt(),  mt,  dk,  yk,  ht,  mot,  dayt) 

DECLARE  SUB  julian . day .numbar  (mt,  dk,  yk,  utt,  jdt(),  djt,  tjt) 

DECLARE  SUB  kaplar  (mt,  act,  aat) 

DECLARE  SUB  output4  (kX,  yrk,  mot,  dak,  tima . of .dayt ,  ltt,  lnt,  htt,  alt,  _ 
azt,  rngt,  langlat,  haadt,  sta.idt,  sat.idt) 

DECLARE  SUB  pos. update  (IX,  elemt(),  timat,  sat.xtO,  sat.xdt(),  elemuf(),  _ 
rt,  aat) 

DECLARE  SUB  rise. sat  (zvi(),  sta.Rt(),  vat,  tOt,  eat,  alemutO) 

DECLARE  SUB  sat. head  (rt,  sat.xt(),  sat.xdtO,  haadt) 

DECLARE  SUB  sat . sphere . coord  (sat.xtO,  rt,  tjt,  tmt,  sat.dct,  sat. lnt,  sat. rat) 
DECLARE  SUB  sta.p. coord  (ltt,  lnt,  htt,  xt()) 

DECLARE  SUB  sta.sat  (sta.xtQ,  sat.xtO,  rngt,  azt,  alt,  langlat) 

DECLARE  SUB  time. update  (kX,  daltatt,  tmt,  jdt,  jdt(),  tjt,  time .of .dayt ,  _ 
mk,  dk,  yk,  ht,  mot,  dayt) 

DECLARE  SUB  yrmoday  (vt,  yrk,  mok,  dyk) 

'COLOR  7,  1 
’SCREEN  0,  1 

DEFDBL  A-Z 
OPTION  BASE  1 

DIM  jdepoch(2) ,  jdstart(2),  jdend(2),  jd(2) ,  sat.x(3),  sat.xd(3),  sta.x(3,  4) 

DIM  elem(9),  elamu(9),  pv(3),  qv(3) ,  sta.z(3),  sta.R(3),  jd . charlie(2) 

’  constants 

CONST  pi  =  3 . 14 1592 653589 79 3t ,  rad  =  pi  /  180t,  twopi  s  pi  ♦  pi 
CONST  j2  =  . 0010826318t 

CONST  nmi.ft  =  ,3048t  /  1852t  ’  n. mi/foot 

CONST  ae  -  6378. 14t  *  earth’s  equatorial  radius  (km.) 

CONST  GE  =  398600. 5t  ’  gao.  grav .  const.  (km)*3/(sec) ‘2 

k  *  60t  *  SQR(GE  /  aa)  /  aa’ (eru)'(3/2)/*in 
Harg.min  =  k 

CONST  fl  =  it  /  298.257 t  '  earth's  flattening  factor 

CONST  a.sqr  -  (2t  -  fl)  *  fl  ’  earth’s  ellipticity  squared 

CONST  we  *  1 . 002737909350795t  *  tsopi  /  1440t’aarth’s  rotation  (radians/minute) 
CONST  sixty  *  60t 

'  Sample  input  (NORAD  Project  Spacatrack  Rpt.  t3,  Dec.  1980): 
amajort  «  "1.04050189":  aect  ■  "0.0086731":  ainclt  «  "72.8435" 
parit  a  "52.6988":  anodet  ■  "115.9689":  manomalyt  *  "110.5714" 
motiont  *  "16.05824518":  dacayt  ■  "0" 
epochdatet  *  "1980.1001":  epochtimet  >  "23.41241138" 
startdatet  *  "1980.1001":  starttimat  ■  "23.41241138" 
anddatat  ■  "1980.1003":  endtimet  »  "0000" 
daltatt  »  "4" 

’  Sample  input  of  arbitrary  higher  altitude,  higher  acc  orbit 
acct  «  "0.0200000":  ainclt  *  "64.94868" 

parit  «  "202.277268":  anodet  «  "55.552464":  manomalyt  «  "126.670896" 
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motion!  *  "12.00000":  decay!  «  "4.281098d-3" 
epochdate!  *  "1983.0201":  epochtime!  «  "00.00" 
startdata!  *  "1983.0201":  starttime!  *  "00.00" 
enddate!  -  "1983.0208":  endtime!  =  "00.00" 
deltat!  *  "IS" 

sta.lat!  *  "3S.12":  eta. Ion!  «  "-8S.12":  at a. alt!  »  "600.0" 

'  Sampla  input  (NAVSPASUR:  Sat.  22222  data  from  SHOVALL,  privata  comm. 
'  Hr.  Frad  Lipp,  July  1988): 
eat. id!  =  "22222" 

acc!  a  "0.0023064":  aincl!  a  "64.94868" 

pari!  a  "202.277268":  anoda!  a  "SS.SS2464":  manoaaly!  «  "126.670896" 

motion!  *  "16.06302":  dacay!  a  "4.281119d-3" 

apochdata!  a  "1983.0201":  apochtima!  a  "00.00" 

startdata!  a  "1983.0201":  starttima!  a  "00.00“ 

anddata!  a  "1983.0203":  andtima!  a  "00.00" 

daltat!  =  "01" 

ata.id!  =  "001":  sta.nama!  a  "Daisy,  Tann." 

sta.lat!  *  "36.12":  ata.lon!  ■  "-8S.12":  sta.alt!  «  "S00.0" 

’  Sample  input  (NAVSPASUR:  Sat.  11111  data  from  SHOVALL,  private  comm. 
’  Mr.  Fred  Lipp,  July  1988): 
sat. id!  =  "11111" 

ecc!  a  "0 . 000S545" :  aincl!  a  "65. 06057" 

peri!  =  "295.41470":  anoda!  a  "272.43497":  manomaly!  a  "258.10682" 

motion!  =  "15.44194":  decay!  a  "0" 

epochdate!  =  "1983.0201":  epochtime!  a  "00.00" 

startdata!  =  "1983.0201":  starttime!  =  "00.00" 

enddate!  =  "1983.0202":  endtime!  =  "00.00" 

deltat!  =  "10" 

sta.id!  =  "001":  sta.name!  =  "Daisy,  Tenn." 

sta.lat!  =  "35.12":  sta.lon!  a  "-8S.12":  sta.alt!  a  "500.0” 

start: 

CLS 

PRINT  "Select  data  input  mode  1,  2,  3  or  4:" 

PRINT 

PRINT  SPC(5) ;  "Keyboard  input  for  satellite  k  ground  station:" 

PRINT  SPC(IO);  "1.  Epoch  of  date." 

PRINT  SPC(IO) ;  "2.  NAVSPASUR  Reference  Data." 

PRINT 

PRINT  SPC(5);  "Data  file  input  for  satellite  k  ground  station:" 

PRINT  SPC(IO);  "3.  Epoch  of  date." 

PRINT  SPC(IO)  ;  "4.  NAVSPASUR  One-Line  Charlie." 
data!  a  "" 

WHILE  data!  <  "1"  OR  data!  >  "4" 
data!  =  INPUT!(1) 

VEND 

CLS 

PRINT  "Select  run  mode  1  or  2:" 

PRINT 

PRINT  SPC(5);  "1.  Create  an  aphemeris." 

PRINT  SPC(5);  "2.  Find  times  satellite  is  above  the  horizon." 

case!  *  "" 

WHILE  case!  <>  "1"  AND  case!  <>  "2" 
case!  *  INPUTS(l) 

VEND 

GOSUB  input .s im. times 

SELECT  CASE  data! 

CASE  "1",  "2" 

CLS 

PRINT  "Satellite  Data:" 

PRINT 
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PRINT  "Epoch  data  (yyyy.mmdd)  7  epochdatet 

LOCATE  3,  29:  INPUT  vt:  IF  vt  *  THEN  vt  »  apochdata$ 
epochdatet  =  vt :  CALL  yrmoday(vt,  yrk,  boA ,  dyk) 

PRINT  "Epoch  time  (hh.mmss)  ?  epochtimet 

LOCATE  4,  29:  INPUT  vt:  IF  vt  «  ""  THE I  vt  •  epochtimet 
epochtimet  *  vt :  epochtime  *  decimal# (vt) 

CALL  julian.day.number(mok,  dyk,  yrk,  epochtime,  jdapochO,  apoch,  tjepoch) 

LOCATE  6,  1:  PRINT  "Input  naan  alaaanta  of  apoch:" 

PRINT  "Satallita  ID  Number  ?  ";  sat. id! 

LOCATE  7,  29:  INPUT  vt:  IF  vt  ■  THEN  vt  *  sat.idt 
sat.idt  =  vt 

PRINT  "Eccantricity  7  acct 

LOCATE  8,  29:  INPUT  vt :  IF  vt  «  ""  THEN  vt  ■  acct 
acct  *  vt:  acc  *  VAL(vt) 

PRINT  "Long,  ascending  noda  (dag)  ?  anodat 
LOCATE  9,  29:  INPUT  vt :  IF  vt  =  ""  THEN  vt  =  anodat 
anodat  =  vt:  anode  -  VAL(vt)  *  rad 

PRINT  "Inclination  (deg)  ?  ";  ainclt 

LOCATE  10,  29:  INPUT  vt:  IF  vt  =  ""  THEN  vt  *  ainclt 
ainclt  =  vt:  aincl  =  VAL(vt)  *  rad 

PRINT  "Arg.  of  perigee  (deg)  ?  parit 

LOCATE  11,  29:  INPUT  vt:  IF  vt  =  ""  THEN  vt  *  perit 
perit  =  vt:  peri  =  VAL(vt)  *  rad 

PRINT  "Mean  anomaly  (deg)  ?  ";  manomalyt 

LOCATE  12,  29:  INPUT  vt :  IF  vt  «  ""  THEN  vt  «  manomalyt 
manomalyt  *  vt:  manomaly  *  VAL(vt)  *  rad 

PRINT  "Mean  motion  (rev/day)  ?  ";  aotiont 

LOCATE  13,  29:  INPUT  vt:  IF  vt  *  ""  THEN  vt  «  motiont 

motiont  =  vt:  motion  *  VAL(vt)  *  twopi  /  1440#’  convert  to  radians/minute 

PRINT  "Decay  (rev/day*2)  ?  decayt 

LOCATE  14,  29:  INPUT  vt:  IF  vt  ■  ""  THEN  vt  *  decayt 
’  convert  to  radians/minute‘2 : 

decayt  =  vt:  decay  =  VAL(vt)  *  twopi  /  (1440#)  *  2 

LOCATE  16,  1:  PRINT  "Press  any  key  to  continue.” 

FOR  iX  =  1  TO  10:  ct  =  INKEYt :  NEXT 
ct  *  INPUTt(l) 

CLS 

LOCATE  1,  1:  PRINT  "Station  Coordinates:" 

PRINT 

PRINT  "Station  Identification  Number  7  " ;  sta.idt 

LOCATE  3,  42:  INPUT  vt:  IF  vt  »  ""  THEN  vt  ■  sta.idt 
sta.idt  s  vt 

PRINT  "Station  Location  or  Name  7  sta.namet 

LOCATE  4,  42:  INPUT  vt:  IF  vt  ■  ""  THEN  vt  ■  sta.namet 
sta.namet  =  vt:  sta.lat  *  VAL(vt)  *  rad 

PRINT  "Station  latitude  (dd.mmss.  South  minus)  7  sta.latt 
LOCATE  5,  42:  INPUT  vt:  IF  vt  «  ""  THEN  vt  ■  sta.latt 
sta.latt  •  vt:  sta.lat  ■  VAL(vt)  *  rad 

PRINT  "Station  longitude  (ddd.mmss.  Vest  minus)  7  sta.lont 
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LOCATE  6,  42:  INPUT  v):  IF  v)  »  ""  THEM  v)  ■  ata.lon) 
ata.lon)  3  v$:  ata.lon  *  VAL(v))  *  rad 

PRINT  "Station  altituda  (feet)  ?  ata.alt) 

LOCATE  7,  42:  INPUT  v):  IF  v)  *  ""  THEM  v)  «  ata.altt 
ata.alt)  *  v):  ata.alt  *  VAL(vl) 

LOCATE  9,  1:  PRINT  "Praas  any  key  to  continoa." 

FOR  iX  3  1  TO  10:  c)  «  INKEYS:  NEXT 
cS  3  INPUTS(l) 

OPEN  "output.dat"  FOR  APPEND  AS  #3 
GOSUB  computa 
CASE  "3" 

OPEN  "atation.dat"  FOR  INPUT  AS  #1 
OPEN  "output.dat"  FOR  APPEND  AS  t3 

DO  UNTIL  EOF(l) 

INPUT  *1,  ata.idS,  ata.latS,  ata.lon),  ata.altS,  ata.name) 
ata.lat  =  VAL(ata.latS)  *  rad 
ata.lon  =  VAL(ata.lon))  *  rad 
ata.alt  3  VAL(ata.altS) 

OPEN  "elementa.dat"  FOR  INPUT  AS  #2 

DO  UNTIL  E0F(2) 

INPUT  t2,  aat.idS,  epochdataS,  epochtime),  accS,  anodaS,  ainclS ,  _ 
peri),  manomaly),  motion),  decay) 

CALL  yrmoday(epochdata),  yrk ,  mok,  dyk) 
epochtime  3  decimalk(epochtime)) 

CALL  julian.day .numbarCnok,  dyk,  yrk,  apochtima,  jdapochO,  epoch, 
tj epoch) 
ecc  3  VAL(acc)) 
anode  3  VAL(anode))  *  rad 
aincl  3  VAL(aincl))  *  rad 
peri  3  VAL(pari))  *  rad 
manomaly  3  VAL(manomaly))  *  rad 

motion  3  VAL(motion))  *  tvopi  /  1440>’  convert  to  radiana/minute 
decay  3  VAL(decayS)  *  twopi  /  (1440k)  *  2  ’  convert  to  rad/min'2 

jd(l)  3  jdatart(l) 
jd(2)  3  jdatart(2) 
jd  3  jd(l)  ♦  jd(2) 
time. of. day  3  atarttima 

GOSUB  computa 

LOOP 
CLOSE  2 
LOOP 

CASE  "4" 

OPEN  "atation.dat"  FOR  INPUT  AS  kl 
OPEN  "output.dat"  FOR  APPEND  AS  k3 

DO  UNTIL  E0F(1) 

INPUT  kl,  ata.id),  ata.lat),  ata.lon),  ata.alt),  ata.name) 
ata.lat  3  VAL(ata.lat))  *  rad 
ata.lon  3  VAL( ata.lon))  *  rad 
ata.alt  3  VAL(eta.alt)) 

OPEN  "elenenta.nea"  FOR  INPUT  AS  k2 
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DO  UNTIL  E0F(2) 

INPUT  82,  ct 

sat. id)  s  MID)(c),  1,  S) 
manomaly)  »  HID)(c),  6,  9) 
motion)  »  HID)(c),  14,  9) 
decay)  »  HID)(c),  22,  7) 

•cc$  «  MID)(c),  28,  9) 
peri)  *  NID)(c) ,  36,  9) 
anode)  *  NID)(e) ,  44,  9) 
aincl)  *  MID)(c) ,  52,  9) 

epochdate)  *  "19"  ♦  NID)(e),  60,  2)  ♦  ♦  NID)(e),  62,  4) 

epochtime)  ■  "00.00" 

CALL  yrmoday(epochdate),  yr* ,  mot,  dyk) 
apochtima  *  decimalt(epochtime)) 

CALL  julian. day .number (mot ,  dyk,  yrk,  apochtima,  jdapochQ ,  apoch,  _ 
tjapoch) 
acc  ■  VAL(ecc)) 
anoda  *  VAL(anode))  *  twopi 
aincl  a  VAL(aincl))  «  tvopi 
pari  -  VAL(peri))  *  tuopi 
manomaly  =  VAL(manomaly$)  *  tuopi 

motion  =  VAL(motion))  *  Herg.min  ’  convart  to  radians /mi nut a 

decay  =  VAL(decay)) 

IF  decay  >  .5#  THEN  decay  *  decay  -  It 

decay  a  decay  *  .00001  *  Harg.min  *  2  ’  convart  to  rad/min*2 

jd(l)  =  jdstart(l) 
jd(2)  =  jdstart(2) 
jd  a  jd(l)  ♦  jd(2) 
time. of. day  a  starttime 

GOSUB  compute 

LOOP 
CLOSE  2 
LOOP 

END  SELECT 
CLOSE 

END 

'  *  *  •  *  *  *  * 

compute: 

'  initialize  times. 

tj  »  tjstart  ’julian  centuries  from  J2000.0 

begin. time  ■  (startdate  -  epoch)  *  1440#  ’minutes 
and. time  *  (enddate  -  apoch)  •  1440#  ’minutes 

elem(l)  »  amajor  ’semi-major  axis 

elem(2)  *  acc  'eccentricity 

elem(3)  *  tzaro  ’time  of  perifocal  passage 

elem(4)  »  aincl  ’inclination  (radians) 

elem(S)  *  anode  ’longitude  of  ascending  node  (radians) 

elem(6)  *  peri  ’argument  of  perifocus  (radians) 

elem(7)  »  motion  ’mean  motion  (revolutions/day) 

elem(8)  ■  manomaly  ’mean  anomaly  (radians) 

elem(9)  ■  decay  'decay  constant  (radians/minute‘2) 

SELECT  CASE  data) 

CASE  "1",  "3" 

’  Since  the  earth  is  treated  as  non-rotating,  elem(5)  must  be  corrected  for 
’  the  accrued  rotation  from  epoch  to  J2000.0 
elem(5)  *  elem(5)  -  158  *  gmst(tjepoch,  08)  *  rad 
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CASE  "2" ,  "4" 

’  NAVSPASUR  one-line  Charlie  correction 

CALL  julian.day .number(l,  1,  1985,  0,  jd. Charlie O ,  dummy,  tj. Charlie) 
wetc  *  (tj epoch  -  tj. Charlie)  *  365268  -  (It  -  gmst(tj .Charlie ,  Ot)  /  248)  _ 
*  .99726956648 

wetc  *  dmod(wetc  *  14408  «  we,  twopi) 
elea(5)  *  elem(5)  ♦  wetc 

elem(5)  *  elem(5)  -  158  *  gmst(tj epoch,  08)  «  rad 
END  SELECT 


CALL  sta.p.coord(eta.lat,  eta. Ion,  eta. alt,  eta.xO) 

FOR  iX  »  1  TO  3 

ata.z(iX)  *  eta.x(iX,  2) 
sta.R(iX)  *  eta.x(iX,  1) 

NEXT 

CALL  pow.npdateCl,  elea(),  begin. tine,  eat.xO,  eat.xdO,  elemuO,  r,  ea) 
CALL  poe  .update(2,  eleaC),  begin. tine,  eat.xO,  eat.xdO,  elemuO,  r,  ea) 

CLS 

PRINT  "Satellite:  sat. id#;  "  Station:  eta. id#;  "  Bta. name# 

PRINT 

SELECT  CASE  case# 


CASE  "1" 

tra  =  starttime  'UT  in  houre 

CALL  apo  .peri  .gee(elem()  ,  elemuO) 

CALL  jd.to.date(jd(),  ml,  dk,  yk,  h8,  mo#,  day$) 

PRINT 

PRINT  day#;  dk;  mot;  yk 

CALL  output4(0 ,  yrk,  mok,  dak,  time. of. day,  eat.lt,  eat. In,  eat.ht,  el,  az ,  _ 
rng,  langle,  head,  sta.id#,  sat. id#) 


’  time-loop 

FOR  time  *  begin. time  TO  end. time  STEP  delta. t 

CALL  poe .update(2,  elemO,  time,  eat.xO,  eat.xdO,  elemuO,  r,  ea) 

'compute  spherical  coordinates 

CALL  sat . sphere .coord(eat .x() ,  r,  t  j ,  tm,  eat.dc,  eat. In,  eat.ra) 

'satellite  heading 

CALL  eat.head(r,  eat.xO,  eat.xdO,  head) 

’earth  ground  track 

CALL  ground. track(r,  eat.dc,  sat.lt,  eat.ht) 

'compute  station-satellite  relations 

CALL  sta.sat(sta.x() ,  eat.xO,  rng,  az,  el,  langle) 

rng  ■  rng  •  ae 

CALL  output4(l,  yk,  mk,  dk,  time. of. day,  sat.lt,  sat. In,  eat.ht,  el,  az,  rng,  _ 
langle,  head,  sta.id#,  eat. id#) 

'update  time 

CALL  time.update(0,  delta. t. hr,  tm,  jd,  jd(),  tj,  time. of. day,  mk,  dk,  yk,  _ 
h8,  mot,  day#) 


NEXT  time 
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PRINT  :  PRINT  "Press  any  key  to  continue." 
FOR  iX  «  1  TO  10:  c$  *  INKEYS :  NEXT 
c$  =  INPUTl(l) 


CASE  "2" 

CALL  output4(0,  yA,  m k,  dA ,  time. of. day,  snt.lt,  set. In,  snt.ht,  el,  nz,  mg 
lnngle,  hend,  sta.idl,  sat. id!) 

tm  =  0 

period  s  twopi  /  motion 

en  *  0 

time  »  begin. time 
initialize: 

CALL  pos  .update(2,  elem(),  time,  sat.xQ,  sat.xdO,  elemu(),  r,  ea) 

CALL  culmnate(sta.zC) ,  we,  time,  ea,  coseta,  elemuO) 

IF  coseta  <  0  THEN 
ea  *  ea  ♦  pi 

CALL  culmnate(sta.z() ,  we,  time,  ea,  coseta,  elemuO) 

END  IF 

time  =  (ea  -  elem(2)  *  SIN(ea))  /  motion  ♦  elemu(3) 
dt.hr  -  (time  -  begin. time)  /  sixty 

CALL  time .update(0 ,  dt.hr,  tm,  jd,  jd(),  t j ,  time. of. day,  mA ,  dA ,  yk,  h# ,  . 
mol,  day!) 

IF  time  >  begin. time  ♦  period  THEN 
time  =  time  -  period 

CALL  time .updated ,  -period,  tm,  jd,  jd(),  t j ,  time. of. day,  mA,  dA,  yk, 
hi,  mol,  dayl) 

GOTO  initialize 
END  IF 

again: 

FOR  iX  *  1  TO  2 

CALL  pos. update(2,  elem(),  time,  sat.x(),  sat.xdO,  elemuO,  r,  ea) 

CALL  culmnate(sta.z() ,  we,  time,  ea,  coseta,  elemuO) 
dt.min  *  (ea  -  elem(2)  *  SIN(ea))  /  motion  ♦  elemu(3)  -  time 
time  s  time  *  dt.min 

CALL  time. updated,  dt.min,  tm,  jd,  jd(),  t j ,  time. of. day,  mA,  dA,  yA,  _ 
ht,  mol,  dayl) 

NEXT 

time.cul  =  time 

CALL  sta.sat(sta.x() ,  sat.x(),  rng,  az,  el,  langle) 

IF  el  >  0  AND  time  >=  begin. time  THEN 

CALL  pos .update(2,  elem(),  time,  sat.x(),  sat.xdO,  elemuO,  r,  ea) 
ea  *  ea  -  .335 
FOR  iX  ■  1  TO  2 

CALL  rise.set(sta.z() ,  sta.RO,  we,  time,  ea,  elemuO) 
dt.min  «  (ea  -  elem(2)  *  SIN(ea))  /  motion  *  elemu(3)  -  time 
time  *  time  *  dt.min 

CALL  time. updated,  dt.min,  tm,  jd,  jd(),  t j ,  time. of. day,  mA,  dA,  _ 
yA,  hi,  mol,  dayl) 

CALL  pos. update(2,  elemO ,  time,  sat.xO,  sat.xdO,  elemuO,  r,  ea) 

NEXT 

time. rise  *  time 

CALL  sta.sat(sta.x(),  sat.xO,  rng,  az,  el,  langle) 

CALL  sat. sphere. coord(sat. x() ,  r,  tj,  tm,  sat.dc,  sat. In,  sat.ra) 

CALL  ground. track(r,  sat.dc,  sat.lt,  sat.ht) 

CALL  sat.head(r,  sat.xO,  sat.xdO,  head) 

CALL  jd.to.date(jd() ,  mA,  dA.  yA,  hi,  mol,  dayl) 

CALL  output4(l,  yA,  mA,  dA,  time. of. day,  sat.lt,  sat. In,  sat.ht,  el,  _ 
az,  rng  •  ae,  langle,  head,  sta.idl,  sat.idl) 
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time.incr  *  dmod(time . cul  -  time. rise,  delta. t) 
tin*  *  time  ♦  time.incr 

CALL  tin*. updat«(l,  tina.incr,  tm ,  jd,  jd(),  tj,  tin*. of. day,  aft,  dk ,  _ 
yk,  h#,  mo#,  day#) 

CALL  poa.updat«(2,  elem(),  tima,  aat.xO,  aat.xdO,  alamuO.  r,  aa) 

CALL  sta.aat(ata.x() ,  aat.x(),  rug,  ax,  al»  langla) 

CALL  sat.  sphere,  coord  (sat.xO,  r,  t j ,  ta,  aat.dc,  aat.ln,  aat.ra) 

CALL  ground. track(r,  aat.de,  aat.lt,  aat.ht) 

CALL  aat.head(r,  aat.xO,  aat.xdO,  head) 

CALL  jd.to.date(jd() ,  mt,  dk,  yk,  hi,  aol,  day#) 

CALL  output4(l,  yk,  ak,  dk,  tiaa.of.day,  aat.lt,  aat.ln,  aat.ht,  al,  _ 
az,  rag  •  aa,  langla,  haad,  at a. id#,  aat.id#) 


DO  WHILE  al  >  0 

tiaa  *  tima  ♦  delta. t 

CALL  time. updated,  delta. t. hr,  ta,  jd,  jd(),  t j ,  tiaa.of.day,  mk,  dk,  _ 
yk,  hS,  ao#,  day#) 

CALL  poa. updated,  alaaQ,  tiaa,  aat.xO,  aat.xdO,  alaau(),  r,  aa) 

CALL  ata.aatCata.xO ,  aat.xO,  rng,  az,  al,  langla) 

IF  al  <  0  THEN 
EXIT  DO 
END  IF 


LOOP 


CALL  sat .sphere. roord(aat . x() ,  r,  t j ,  ta,  sat.dc,  sat. In,  sat.ra) 

CALL  ground. trac.  r,  sat.dc,  sat.lt,  sat.ht) 

CALL  sat.head(r,  sat.xO,  sat.xdO,  haad) 

CALL  jd.to.date(jd() ,  mk,  dk,  yk,  ht,  mo#,  day#) 

CALL  output4(l,  yk,  mk,  dk,  tiaa.of.day,  sat.lt,  sat. In,  sat.ht,  al,  _ 
az,  rng  *  aa,  langla,  head,  sta.id#,  sat. id#) 


save. tima  *  tima 
FOR  i*  =  1  TO  2 

CALL  pos. update (2,  alam(),  time,  aat.xO,  sat.xdO,  alamuO,  r,  aa) 

CALL  risa.sat(sta.z() ,  sta.RO,  we,  tima,  aa,  alamuO) 
dt.min  *  (aa  -  elem(2)  •  SIN(ea))  /  motion  ♦  elemu(3)  -  tima 
tima  *  tima  ♦  dt.min 

CALL  tima .updataCl ,  dt.min,  tm,  jd,  jd(),  t j ,  tima. of. day,  mk,  dk,  yk,  _ 
h«,  mo#,  day#) 

NEXT 

CALL  pos . updated,  elemO,  tima,  aat.xO,  aat.xdO,  alamuO,  r,  aa) 

CALL  sta.sat(sta.x() ,  sat.xO,  rng,  az,  al,  langla) 

CALL  sat . sphere . coord (sat ,x() ,  r,  t j ,  tm,  sat.dc,  aat.ln,  sat.ra) 

CALL  ground. track(r,  sat.dc,  sat.lt,  sat.ht) 

CALL  sat.headCr,  sat.xO,  sat.xdO,  haad) 

CALL  jd.to.date(jd() .  mk,  dk,  yk,  hi,  mo#,  day#) 

CALL  output4(l,  yk,  mk,  dk,  tima. of. day,  sat.lt,  sat. In,  sat.ht,  al,  . 
az,  rag  *  aa,  langla,  haad,  sta.id#,  aat.id#) 

PRINT 


END  IF 

tima  >  tima. cul  *  period 

dt.hr  ■  (tima  -  begin. tima)  /  aixty  -  tm 

CALL  time. updated,  dt.hr,  ta,  jd,  jd(),  t j ,  time. of. day,  mk,  dk,  yk,  ht,  _ 
mo#,  day#) 

IF  tima  <  and. tima  THEN  GOTO  again 


END  SELECT 
RETURN 

input. aim. timaa: 

CLS 

PRINT  "Simulation-Time  Data:" 

PRINT 

PRINT  "Starting  data  (yyyy.mmdd)  ?  ";  atartdate# 
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LOCATE  3,  29:  INPUT  v$:  IF  vl  »  1111  THEN  vl  •  startdatal 
startdatal  «  v$:  CALL  yraoday(vl,  yrk,  mok,  dyk) 

PRINT  "Starting  tin*  (hh.mmss)  ?  start time* 

LOCATE  4,  29:  INPUT  vl:  IF  vl  *  ""  THEN  vl  ■  starttimal 
starttimal  *  vl:  *t art tin*  ■  d*cimal#(vl) 

CALL  j alias. day . number (aok,  dyk,  yrk,  starttima,  jdatartO,  atartdat*,  _ 

tjatart) 

jd(l)  =  jdstart(l) 
jd(2)  =  jdstart(2) 

jd  *  jdCl)  ♦  jd(2) 
tin*. of .day  *  atarttim* 

PRINT  "Ending  dat*  (yyyy.mmdd)  ?  anddatal 

LOCATE  5,  29:  INPUT  vl:  IF  vl  *  ""  THEN  vl  *  anddatal 
anddatal  =  vl:  CALL  yrmoday(v|,  yrk,  aok,  dyk) 

PRINT  "Ending  tim*  (hh.mmaa)  ?  " ;  *ndtim*$ 

LOCATE  6,  29:  INPUT  vl:  IF  vl  ■  ""  THEN  vl  •  *ndtim*l 
•ndtimel  =  vl:  sndtim*  *  d*cinalt(v|) 

CALL  julian.d&y .numb*r(mok,  dyk,  yrk,  andtia*,  jd*nd(),  enddata,  tjand) 

PRINT  "Tima  incrament  (minutaa)  ?  daltatl 

LOCATE  7,  29:  INPUT  vl :  IF  vl  *  ""  THEN  vl  »  daltatl 

daltatl  =  vl:  dalta.t  =  VAL(vl):  dalta.t.hr  =  dalta.t  /  sixty 

LOCATE  9,  1:  PRINT  "Press  any  key  to  continua." 

FOR  if.  =  1  TO  10-  .t  -  INKEYI :  NEXT 
cl  =  INPUTl(l) 

RETURN 

>*****«**«*«•*»**•*«*«****«***• ***************************************** **»**•» 
FUNCTION  acos#  (x#)  *  02-10-88  Rav.  02-24-88 

’  Th*  arccosin*  function  darivad  from  th*  ATN  function  aith  no  singularitias . 

'  The  angl*  is  raturnad  in  th*  (0,pi)  intarval. 

’Not*  that  in  ralational  comparisons,  -1  ia  "trua"  and  0  is  "falsa". 

CONST  pi  =  3.141592653589793#,  aps  a  ID-33 
IF  ABS(xt)  <=  1*  THEN 

acos#  =  ATN(SQR( 1#  -  x#  *  x#)  /  (x#  -  aps  *  (x#  =  0#)))  -  pi  *  (x*  <  0#) 

ELSE 

PRINT 

PRINT  "Error  in  acos  function.  ABS(arg)  >  1" 

PRINT  "arg  =  ";  x# 

PRINT 
END  IF 

END  FUNCTION 


SUB  apo.pari.ga*  (*l*m#(),  *l*mu#())  *  09-14-88.  Rav  09-14-88. 

’  Comput*  location  and  haight  of  apogaa  and  pariga*  abov*  oblat*  aarth. 

’  NOTE:  (1)  Th*  output  from  this  routin*  ia  approximataly  valid  for  th* 

'  currant  orbit  only.  In  particular,  th*  longitud*  is  not  corractad 

’  for  th*  aarth'a  rotation  during  th*  currant  orbit.  Also,  th* 

’  ascending  nod*  and  argumant  of  pariga*  ara  not  corractad  for 

’  thair  rat as  of  chang*  during  th*  currant  orbit. 

’  (2)  For  naarly  circular  orbits,  tha  haight  of  th*  satallit*  sill 

’  probably  gat  losar  than  th*  pariga*  valu*  and  highar  than  th* 

’  apogaa  valu*  at  locations  other  than  th*  pariga*  and  apogaa 

’  positions— this  is  an  affect  of  th*  earth’s  oblatanass. 
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CONST  pi  *  3 . 141592653589793*,  rad  «  pi  /  180#,  tsopi  *  pi  ♦  pi 


'  spherical  sarth  values: 

apogee#  *  elemu#(l)  *  (1#  -  elen#(2)) 

perigee#  »  elemu#(l)  *  (1#  ♦  elem#(2)) 

'  longitudes  and  declinations: 
sincl#  «  SIN(elem#(4)) 
speri#  «  SIN(eleau#(6)) 
apo.dc#  a  asin#(speri#  *  sincl#) 

apo.ln#  *  eleau#(5)  +  atan2#(C0S(elea#(4))  •  speri#,  C0S(eleau#(6) )) 
apo.ln#  a  daod* (apo.ln#,  tsopi) 

IF  apo.ln#  >  pi  THEN  apo.ln#  ■  apo.ln#  -  tsopi 

’  next  are  approxiaations  should  update  by  peri-  k  nodal-rates 
'  for  1/2  revolution: 
peri.dc#  ■  -apo.dc# 
peri. In#  *  daod Capo . In#  ♦  pi,  tsopi) 

IF  peri. In#  >  pi  THEN  peri. In#  ■  peri. In#  -  tsopi 

CALL  ground. track(apogee# ,  apo.dc# ,  apo.lt#,  apo.ht#) 

CALL  ground. trackCperigee#,  peri.dc#,  peri. It#,  peri.ht#) 

PRINT 

PRINT  "Apogee:" 

PRINT  USING  "*  *#*#.#  "Height  a  apo.ht#;  "  km.  " 

PRINT  USING  "#  ###.##*#";  "Latitude  =  apo.lt#  /  rad 

PRINT  USING  "t  ##«.#«<#";  "Longitude  a  daodCapo.ln#  /  rad,  360#) 

PRINT 

PRINT  "Perigee:" 

PRINT  USING  "k  ##«#.#  k" ;  "Height  =  peri.ht#;  "  km.  " 

PRINT  USING  "k  ###.#«#";  "Latitude  «  peri. It#  /  rad 

PRINT  USING  "*  ###.####";  "Longitude  a  dmod(peri.ln#  /  rad,  360#) 

PRINT 

IF  apo.ht#  <  0  THEN 

PRINT  "WARNING:  Satellite  is  suborbital  and  sill  impact  earth." 

PRINT  "Do  you  sant  to  continue?  Yes  or  No." 

PRINT 

WHILE  c$  <>  "y"  AND  c#  <>  "n" 
c$  *  LCASE8 ( INPUTS ( 1 ) ) 

WEND 

IF  c$  a  "n"  THEN  END 
END  IF 
END  SUB 

>•**••••**•*•«***•*••*•*•*•*»**•*•****•**•******«*•*•**•••*••*******•******•**• 
FUNCTION  asin#  (x#)  ’  02-10-88  Rev.  02-24-88 

»  The  arcsine  function  derived  from  the  ATN  function  sith  no  singularities. 

'  The  angle  is  returned  in  the  (-pi, pi)  interval. 

'Note  that  in  relational  comparisons,  -1  is  "true"  and  0  is  "false". 

CONST  eps  a  ID-33 
IF  ABS(x#)  <■  1#  THEN 

aein#  »  ATN(x#  /  (SQR(1#  -  x#  *  x#)  -  eps  *  (ABS(x#)  ■  1#))) 

ELSE 

PRINT 

PRINT  "Error  in  asin  function.  ABS(arg)  >  1" 

PRINT  "arg  «  ";  x# 

PRINT 
END  IF 
END  FUNCTION 


FUNCTION  atan2#  (y# ,  x#)  •  02-10-88  Rev.  02-24-88 
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’  The  two  argument  quadrant  determining  arctangent  function. 

1  The  angle  is  returned  in  the  (-pi, pi)  interval. 

’Note  that  in  relational  comparisons,  -1  is  "true"  and  0  is  "false". 
CONST  pi  -  3.141592663589793#,  eps  =  ID-33 

atan2#  *  ATN(y#  /  (x#  -  eps  *  (x#  »  0*)))  -  pi  *  (x#  <  0#)  *  (SGN(y#)  _ 
-  (y#  *  0#)) 


END  FUNCTION 


SUB  culmnate  (zvi(),  net,  tOt,  eat,  cosetat,  elemutQ) 


Compute  eccentric  anomaly  for  the  time  of  culmination  of  a  satellite 
01-30-89.  Rev.  02-02-89  C  1000. 


NOTE:  This  particular  form  of  the  culminate  function  was  written  from  the 

development  of  30  Jain  89.  It  assumes  a  non-rotating  earth  so  that  the 
station  longitude  is  fixed.  Rotation  is  accomplished  by  subtracting  the 
accumulated  rotation  from  the  right  ascension  of  the  ascending  node  of 
the  satellite.  This  is  the  method  used  to  compute  the  ephemeris  and  is 
also  equivalent  to  the  method  used  by  NAVSPASUR. 

lambda  =  station  geodetic  longitude 

phi  =  station  geodetic  latitude 

we  -  earth’s  rotation  (sidereal  rate  of  change) 

tO  =  time  the  elements,  elemuO,  sere  last  updated 

capt  -  T  —  time  of  latest  perifocal  passage 

mot  =  mean  motion 

ea  =  eccentric  anomaly  E 

ecc  =  eccentricity  e 

amjr  =  semimajor  axis  (cancels  out  in  all  formulas,  so  not  included) 


CONST  pi  *  3.141592653589793#,  rad  *  pi  /  180#,  twopi  »  pi  ♦  pi 


ecc#  =  elemu#(2) 
capt*  *  elemu#(3) 
incl#  =  elemu#(4) 
node#  »  elemu#(5) 
peri#  *  elemu#(6) 
mot#  =  elemu#(7) 


DIM  pv#(3) ,  pvd#(3) ,  pvdd#(3) ,  qv#(3),  qvd#(3),  qvdd#(3) 

DIM  rvec*(3),  rvecd#(3) ,  rvecdd#(3) 

’  zv  =  Z-vector,  unit  normal  at  the  station 

’  pv  =  P-vector,  qv  *  Q-vector  (first  two  columns  of  the  Euler  matrix) 

’  rvec  =  R-vector  (satellite  position),  rvecd  ■  dR/dE,  rvecdd  =  d2R/dE2 

’  rv  *  r - magnitude  of  the  R-vector,  rvd  ■  dr/dE,  rvdd  «  d2r/dE2 

’  gamma  ■  node  -  accrued  rotational  motion 

cosn#  *  C0S(peri#) 
sinv#  s  SIN(peri#) 
cosi#  a  C0S(incl#) 
sini#  *  SIN(incl#) 


corr#  «  1 

WHILE  ABS(eorr#)  >  .0001 

cosea#  >  C0S(ea#) 
ecosea#  ■  ecc#  *  cosea# 
sinea#  «  SIN(ea#) 
esinea#  =  ecc#  •  sinea# 

gamma#  a  node#  -  «e#  •  ((ea#  -  esinea#)  /  mot#  *  capt#  -  to#) 
gammad#  «  -ve#  *  (1*  -  ecosea#)  /  mot# 
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gammadd#  «  -Be#  *  eainea#  /  mot# 

cosg#  =  C0S(gamma#) 

•  ingf  =  SIN(gamma#) 

rtaccf  »  SQR( l#  -  ecc#  *  acct) 

r»#  ■  1#  -  ecoaea# 
rvd#  «  eainea# 
rvddt  «  acosaa# 


xv#  «  coaaa#  -  acct 
XBdS  a  -ainea# 
xsddt  a  -coaaa# 
ya #  a  rtacc#  «  ainaa# 
yed#  a  rtacc#  *  coaaa# 
ysdd#  *  -yw# 

pv#(l)  a  coav#  *  coag#  -  ainB*  *  aing#  *  coai# 
pv#(2)  a  Cobb*  *  aing#  ♦  aine#  *  coag#  a  coai# 
pv#(3)  =  Bins#  *  aim# 

qv#(l)  a  -ainB#  *  coag#  -  coau#  *  aing#  a  coai# 
qv#(2)  a  -ains#  *  aing#  ♦  cobb#  a  coag#  a  coai# 
qv#(3)  a  cobb#  a  aini# 

pvd#(l)  a  -pv#(2)  a  gammad# 
pvd*(2)  =  pv#(l)  a  gamaad# 
pvd#(3)  =  0# 

qvd#(l)  =  -qv#(2)  a  gammad# 
qvd#(2)  =  qv*(l)  *  gammad# 
pvd#(3)  =  0# 


pvdd#(l)  =  -pvd*(2)  *  gammad#  -  pv#(2)  a  gammadd# 
pvdd#(2)  =  pvd#(l)  a  gammad#  ♦  pv#(2)  a  gammadd# 
pvdd#(3)  *  0# 

qvdd#( 1)  a  -qvd#(2)  »  gammad#  -  qv#(2)  a  gammadd# 
qvdd#(2)  *  qvd#(l)  a  gammad#  +  qv#(2)  a  gammadd# 
pvdd*(3)  =  0# 


FOR  i%  a  i  TO  3 

rvec#(iX)  *  xb#  a  pv#(i%)  ♦  ya#  a  qV#(iX) 

rvecd#(iX)  *  XBd#  a  pv#(iX)  ♦  xb#  a  pvd#(iX)  ♦  yvd#  *  qv#(iX) 

♦  yB#  a  qvd*(iX) 

temp#  a  xsdd#  a  pv#(iX)  ♦  2#  a  Xad#  a  pvd#(iX)  ♦  xb#  a  pvdd#(iX) 
rvecdd#(iX)  =  temp#  ♦  yadd#  a  qv#(iX)  ♦  2#  *  yad#  a  qvd#(iX) 

♦  ya*  a  qvdd#(iX) 

NEXT 


coaeta#  *  dot#(zv#(),  rvec#0)  /  rv# 

fe#  a  dot#(zv#( ) ,  rvecd#())  -  rvd#  a  coaeta# 

rrd#  »  rvd#  /  rv# 

fed#  a  dot#(zv#() ,  rvecdd#())  ♦  coaeta#  *  (rrd#  -  rvdd*) 
fed#  ■  fed*  -  rrd#  a  dot#(zv#(),  rvecd#()) 

corr#  a  fe#  /  fed# 

IF  corr#  >  pi  THE* 

ea#  ■  ea#  ♦  tsopi 

ELSE 

ea*  »  ea#  -  corr# 

IF  coaeta#  <  0  THEN  ea#  a  «a#  ♦  pj 
END  IF 


VEND 
END  SUB 
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* ******* ********************************************************** ************* 
FUNCTION  decimal#  (v$)  »  03-21-88.  Rev  03-21-88. 

'  Convert  ddd. massif  or  hh.mmssfff  to  decimal  degrees  or  decimal  hours. 

ix  =  INSTR(v$,  ".") 

IF  ix  a  o  THEN 

dscimalt  *  VAL(v$) 

ELSE 

x#  *  VAL(LEFT3(v$,  ix)) 
sn  a  0 

IF  x#  <  0#  THEN 
sn  s  -1 

X*  *  -7 

END  IF 

«$  «  v$  ♦  "0000" 

y#  «  VAL(MID$(s$,  ix  ♦  1 ,  2)) 

z#  *  VAL(MID$(s$,  ix  ♦  3,  2)  ♦  ♦  RIGHT$(b$,  LEN(b$)  -  ix  -  4)) 

xt  «  (x#  /  60#  ♦  y#)  /  60#  ♦  x# 

IF  sn  THEN  x#  =  -x# 
decimal#  =  x# 

END  IF 

END  FUNCTION 

'*•****•******•*•*****••****«•*•***•*••******•**•**••••••**«•**••**•••******•*• 

FUNCTION  Imod#  (x#,  m#)  '  02-10-88  Rev.  02-24-88 

1  Provides  x  MOD  m  lor  real-valued  functions. 

IF  m#  <>  0#  THEN 

dmod#  =  x#  -  m#  *  INT(x#  /  m#) 

ELSE 

PRINT 

PRINT  "Error  in  dmod  function.  Modulus  cannot  be  xero." 

PRINT 
END  IF 

END  FUNCTION 

■a**************************!*************************************************** 

FUNCTION  dot#  (a#(),  b#()) 

'  Dot  product  of  two  vectors. 

dot#  =  a#( 1)  *  b#(l)  ♦  a# (2)  *  b#(2)  ♦  a#(3)  *  b#(3) 

END  FUNCTION 

I****************************************************************************** 

SUB  euler  (peri#,  node#,  incl#,  psn#(),  vel*()) 

'  09-05-88.  Rev.  01-16-89. 

’Compute  the  first  tvo  columns  (p  k  q)  of  the  Euler  matrix. 

'Rotate  in-plane  position  and  velocity  components  to  rectangular  equatorial 
’  coordinates. 

'psnO  and  vel()  are  both  input  and  output  vectors.  It  is  assumed  that  the 
’  third  component  of  both  psn()  and  vel()  is  xero  on  input. 

cat  a  C0S(peri#) 

SB#  a  SIN(peri#) 

cn#  a  C0S(node#) 

sn#  *  SIN(nodeS) 

ci#  *  C0S(incl#) 

si#  a  siNCincl#) 

px#  »  cb*  *  cn#  -  sb#  *  sn#  *  ci# 
py*  «  cb#  *  sn#  ♦  sb#  *  cn*  *  ci# 
px#  ■  SB#  *  Si# 

qx#  *  -sb*  *  cn#  -  cw #  *  sn*  •  ci# 

qy#  «  -sb#  »  sn#  *  cu#  *  cn#  *  ci# 

qx#  ■  cb#  *  si# 
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tl#  *  psn#(l) 
t2#  =  psn#(2) 

psn#(l)  =  px#  *  tl#  ♦  qx#  *  t2# 

psn#(2)  »  py#  *  tl#  ♦  qy#  *  t2# 

psn#(3)  »  px#  *  tl#  *  qz#  *  t2# 

tl#  a  vel#(l) 
t2#  a  vel#(2) 

v»l*(l)  a  px#  *  tl#  ♦  qx#  *  t2# 

Vel#(2)  a  py#  *  tl#  +  qy*  *  t2# 

Vel*(3)  a  pz#  *  tl#  ♦  qz#  *  t2# 

END  SUB 

>  ****************************************************************************** 
FUNCTION  gmst#  (tu#,  tm#)  *  03-03-88.  Rev  03-03-88. 

'  Comput •  Greenwich  mean  sidereal  tin*. 

»  tu#  a  number  of  centuries  of  36525  day*  of  universal  tima  alapaad  since 
»  2000  Jan  1,  12h  UT1  (JD  2461646  UT1). 

•  tm#  a  tima  of  day 

’  Tha  Astronomical  Almanac,  1984,  pp  S13-S15. 

’  Almanac  for  Computers  1988,  pp  B2-B3. 

CONST  pi  =  3.141592653589793#,  rad  a  pi  /  180# 

CONST  cO  =  24110.54841#,  cl  =  8640184.812866#,  c2  =  9 . 310400000000001D-02 
CONST  c3  =  -.0000062#,  hourpersee  =  1#  /  3600#,  meantoapp  =  -.00029# 

CONST  dO  =  125.004452#,  dl  =  -.0529538#,  d2  =  .002071# 

t#  =  ( ( ( c3  *  tu#  +  c2)  *  tu#  +  cl)  *  tu*  ♦  cO)  *  hourpersee  +  tm#’  GMST 
'The  following  two  lines  sill  convert  GMST  to  CAST 
’lunarnode#  =  (d2  *  tu#  ♦  dl)  *  tu#  *  tO# 

*t*  ■  t#  ♦  meantoapp  *  SINClunarnode#  *  rad) 
gmst#  a  dmod# (t# ,  24#) 

END  FUNCTION 

>•**••••*•****»***•*****• ***************************************** ************* 
SUB  ground. track  (r#,  sat.dc#,  sat. It#,  sat.ht#)  '09-14-88.  Rev  09-14-88. 

CONST  ae  =  6378.14#  '  earth’s  equatorial  radius  (km.) 

CONST  fl  =  1*  /  298.257#  '  earth’s  flattening  factor 

CONST  e.sqr  =  (2#  -  fl)  *  fl  ’  earth’s  allipticity  squared 

CONST  12  =  (1*  -  fl)  »  (1#  -  fl) 

rkm#  a  -#  *  ae 
ps#  *  sat.dc# 
al:  cp#  *  C0S(ps#) 

rc#  a  aa  *  SQR((1#  -  e.sqr)  /  (1#  -  e.sqr  *  cp#  *  cp#)) 
rcr#  *  rc#  /  rkm# 

sat. It#  a  ATN(TAN(ps#)  /  12)  '  satellite  latitude 

sn#  >  SINCsat.lt#  -  ps#) 

hr#  »  SQR(1#  -  rcr#  *  rcr#  *  sn#  *  sn#)  -  rcr#  *  C0S(sat.lt#  -  ps#) 
pn#  >  sat.dc#  -  asin(hr#  •  sn#) 

IF  (ABS(pn#  -  ps#)  >  .0000001)  THEN  ps*  >  pn#:  GOTO  el 
sat.ht#  a  hr#  •  rkm#  ’  satellite  height 

END  SUB 


SUB  hr. min. sec  (xS,  hrX,  minX,  sec#,  n%) 

’  Convert  decimal  hours  to:  hh  mm  ss  if  it  *  0, 

'  hh  mm  ss.f  if  nX  a  i, 

’  hh  mm  ss.ff  if  nX  *  2, 

'  hh  mm  ss.fff  if  nX  =  3. 

CONST  to  ■  1#,  tl  »  10#,  t2  a  100#,  t3  ■  1000# 


CONST  rO  *  1#  /  7200# ,  rl  ■  rO  /  tl ,  r2  »  rl  /  t2 .  r3  «  r2  /  t3 


SELECT  CASE  nX 
CASE  0 

round  «  rO 
trunc  *  tO 
CASE  1 

round  *  il 
trunc  ■  tl 
CASE  2 

round  *  r2 
trunc  =  t2 
CASE  3 

round  =  r3 
trunc  *  t3 
CASE  ELSE 
PRINT 

PRINT  "Error  in  hr. min. sec  subroutine.  Must  have  0  <=  nX  <=  3." 

PRINT  "nX  =  nX 

PRINT 

EXIT  SUB 

END  SELECT 

v#  =  ABS(x#)  ♦  round 
hrX  =  INT(vf) 
v#  =  (v#  -  hrX)  *  60# 
minX  =  INT(v#) 

see#  =  INT((60#  *  (v#  -  minX)  *  trunc))  /  trunc 
END  SUB 

>  ****************************************************************************** 
SUB  jd. to. date  (jd#(),  mA,  dA,  yA,  h#,  mol,  dayl)  '  02-24-88  Rev.  03-02-88 
’  Convert  Julian  Day  Number  to  month  number,  day,  year,  hour,  month  name  and 
1  day  of  the  week.  Limited  to  the  Gregorian  Calendar. 

’  Fliegel  A  VanFlandern,  Comm.  ACM,  Vol.  11,  No.  10,  Oct.  1968,  pg.  657. 

CONST  wl  =  "HonTueWedThuFriSatSun" 

CONST  ddl  *  " JanFebMarAprMayJunJulAugSepOctNovDec" 

j4A  =  INT(jd#(l)  *  jd#(2)  ♦  .5) 

h#  =  (( jd#(l)  -  j4A)  +  jd*(2)  ♦  .6)  *  24# 

1A  =  INT(j4A  ♦  68569) 
nA  =  I  NT  (4  *  1A  /  146097) 

1A  »  1A  -  INT(( 146097  *  nA  ♦  3)  /  4) 
yA  =  I NT (4000  *  (1A  ♦  1)  /  1461001) 

1A  *  1A  -  IMT(1461  *  jA  /  4)  ♦  31 
mA  «  INT(80  •  1A  /  2447) 
dA  »  1A  -  IHT(2447  *  mA  /  80) 

1A  *  INTCmA  /  11) 

mA«mA  +  2-12*lA 

yA  »  100  •  (nA  -  49)  ♦  yA  ♦  1A 

vnX  »  j4A  -  7  *  INT(j4A  /  7)  ♦  1 
day!  »  MIDl(vl,  3  *  wnX  -  2,  3) 
mol  a  MIDKddl,  3  *  mA  -  2,  3) 

END  SUB 


SUB  julian. day. number  (mA,  dA,  yA,  ut<,  jd#(),  dj#,  tj#) 

*  02-25-88  Rev.  03-23-88 

’  Conversion  of  calendar  date:  mA  >  month,  dA  *  day,  yA  ■  4-digit  year, 

’  and  ut#  «  Universal  Time,  to  Julian  Day  Number.  Year  must  be  between  1801 
'  and  2099. 
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’  Almanac  for  Computers,  1988,  pg.  B-2. 


Note  that  jd#  is  a  subscripted  double  precision  array.  At  any  time,  the 
Julian  Day  Number  is  given  by  jdt(l)  ♦  jd#(2). 

For  ease  of  programming,  the  user  may  put  the  entire  epoch  in  jd#(l)  and  set 
jdt(2)  =  0#. 

For  maximum  accuracy,  set  jdf(l)  ■  the  most  recent  midnight  at  or  before  the 
epoch  and  set  jd#(2)  *  the  fractional  part  of  a  day  elapsed  between  jd#(l) 
and  the  epoch.  As  an  example,  compute  the  number  of  daye  from  epoch 
J2000.0  *  JD  24S1545.0  using  the  code: 

days#  =  ( jd*(l)  -  2451545.0#)  ♦  jd#(2) 

Note,  particularly,  the  extra  pair  of  parenthesis. 

dj#  =  days  and  fraction  from  J2000.0.  tjt  »  julian  centuries  from  J2000.0 


CONST  cl  =  1721013.5#,  c2  *  190002.5#,  jcent  «  36525# 

IF  (1801  <=  yk)  AND  (yk  <*  2099)  THEN 

jd#(l)  =  367  *  yk  -  INT(7  *  (yk  ♦  INT((mk  ♦  9)  /  12))  /  4)  _ 

+  INT( (275  *  mk)  /  9)  +  dk  ♦  cl  -  .5#  *  SGN(100  *  yk  ♦  mk  -  c2)  +  .5# 
jd#(2)  =  ut#  /  24 

dj#  =  ( jd#( 1)  -  2451545#)  ♦  jd#(2) 
tj#  =  dj#  /  jcent 
ELSE 

PRINT 

PRINT  "Error  in  julian .day .number .  Year  is  not  between  1801  and  2099." 

PRINT  "Year  =  " ;  yk 
PRINT 
END  IF 
END  SUB 

■ I**************.**.*********************************************************** 

SUB  kepler  (m*.  ec* ,  ea#)  ’  02-25-88.  Rev  02-26-88 

’High  precision,  non-iterative,  solution  to  Kepler's  equation  for  the  ellipse. 
’Accuracy  is  10E-18  using  all  terms  or  10E-15  if  the  last  term  is  omitted. 
’Although  not  as  compact  as  the  Nevton-Raphson  procedure,  one  square  root,  one 
’cube  root,  and  only  two  trigonometric  functions  are  evaluated  thus  making  this 
’procedure  extremely  fast.  Can  be  extended  to  the  hyperbolic  case  -  see  Ref. 
'Ref.:  S.  Nikkalo ,  Cel.  Hech.,  Vol.  40,  1987,  pp  329-334. 

’m#  =  mean  anomaly,  ec#  *  eccentricity  and  ea#  *  eccentric  anomaly. 

CONST  pi  =  3.141592653589793#,  tvopi  *  pi  ♦  pi,  third  *  1#  /  3* 

CONST  r3  =  1#  /  6#,  r4  «  1#  /  24#,  r5  *  1#  /  120# 

IF  ec#  >*  1#  OR  ec#  <  0#  THEN 
PRINT 

PRINT  "Error  in  Kepler.  The  eccentricity  is  not  in  the  [0,1)  interval." 
PRINT  "Eccentricity  ■  ec# 

PRINT 

ELSE 

m#  *  m#  -  tvopi  *  INT(m#  /  tvopi)’  assure  that  m#  is 

IF  m#  >  pi  THEN  m#  «  m#  -  tvopi’  in  the  (-pi, pi)  interval 

•Compute  the  initial  approximation: 
g#  «  1  /  (4  *  ec#  ♦  .5) 
at  *  (1  -  ec*)  *  gt 
b#  =  .5  *  m#  *  g# 

arg#  *  ABS(b#)  ♦  SQR(b#  *  b#  +  a#  *  a#  *  at) 
z#  *  arg#  *  third 
IF  b#  <  0  THEN  z#  *  -z# 

s#  »  z#  -  a#  /  z# 

s#  «  s#  -  .078  *  s#  *  5  /  (1  ♦  act) 

ea#  ■  m#  ♦  ec#  *  s#  •  (3  -  4  *  s#  •  s#)’  abs(dea/ea)  <*  0.002 
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’Refine  the  approximation: 
f  2t  »  act  *  SIN(eat) 
f 3#  a  act  a  COS(aaf) 
fO#  a  a at  -  12 t  -  at 
lit  a  It  -  13t 
14t  a  -12t 
15t  a  -13t 

’This  next  correction  tarn  ia  sufficient  il  acc  <  0.2S 
dt  a  -lot  /  lit 

'Use  aora  ol  tha  following  taraa  to  increase  accuracy 
dt  a  -lot  /  (lit  ♦  .St  *  12t  a  dt) 

dt  a  -fot  /  (lit  ♦  dt  a  (.St  •  !2f  ♦  dt  a  r3  a  13t)) 

dt  =  -lot  /  (lit  ♦  dt  a  (.St  a  I2t  ♦  dt  *  (r3  *  13t  ♦  dt  *  r4  *  14t))) 

dt  a  -fot  /  (lit  ♦  dt  a  (.St  a  I2t  ♦  dt  a  (r3  a  I3t  ♦  dt  _ 

a  (r4  a  I4t  Ht  •  rS  a  ISt)))) 
aat  a  aat  a  dt 
END  IF 
END  SUB 

'••••••••♦•••••a*************************************************************** 

SUB  output4  (kX,  yrk ,  mot,  dak,  tina.ol.day,  It,  Is,  ht,  al,  az,  rng,  _ 
langla,  haad,  sta.idl,  sat.idl) 

CONST  pi  a  3 . 141592653589793t ,  rad  «  pi  /  180t,  tp  =  pi  ♦  pi 
c$  =  CHR$(34) 

SELECT  CASE  kX 
CASE  0 

PRINT  "  height  " ; 

PRINT  "  look" 

PRINT  "year  mo  da  hr  an  eec  lat  long  (km)  elev 

PRINT  "azia  range  angle  head" 

PRINT 

CASE  1 

CALL  hr.min.aec(tirne.ol .day#,  hrX,  ainX,  aact,  1) 

It  a  It  /  rad 

In  a  dmodQn  /  rad,  360t) 

la  =  langla 

hd  =  haad  /  rad 

PRINT  USING  "tttt  tt  tt  tt  tt  tt.t  yrk;  mot;  dak;  hrX;  minX;  sec; 
PRINT  USING  "ttt.tt  ttt.tt  ttttt.t  ttt.tt  ";  It;  In;  ht;  el; 

PRINT  USING  "ttt.tt  ttttt  tt.tt  ttt.tt";  az;  rng;  la;  hd 
PRINT  t3 ,  USING  "!k!  !k!  ";  d;  sta.idl;  cl;  c$;  sat.idl;  c$; 

PRINT  t3 ,  USING  "ttttt!  tttt  ";  d;  yrk;  cl;  cl;  aok;  cl; 

PRINT  t3,  USING  "tttt  tttt  ";  c$;  dak;  cl;  el;  hrX;  cl; 

PRINT  t3 ,  USING  "tttt  ttt.tt  cl;  minX;  cl;  cl;  sec;  cl; 

PRINT  t3 ,  USING  "tttt.tt!  tttt.tt!  ";  cl;  It;  cl;  cl;  In;  cl; 

PRINT  t3,  USING  "tttttt.t!  tttt.tt!  cl;  ht;  cl;  cl;  al;  cl; 

PRINT  t3,  USING  "tttt.tt!  tttttt!  ";  cl;  az;  cl;  d;  rng;  cl; 

PRINT  t3,  USING  "ttt.tt!  tttt.tt!";  cl;  la;  cl;  cl;  hd;  cl 
CASE  ELSE 
PRINT 

PRINT  "Error  in  output4  subroutine.  Rust  have  0  <■  kX  <=  1." 

PRINT  "kX  «  ";  kX 
END  SELECT 
END  SUB 


» 


SUB  pos. update  (IX,  eleatQ,  tiaat,  sat.xt() 
STATIC 


sat.xdt(),  alaaut(),  rt,  aat)  _ 
09-07-88.  Rev.  03-09-89. 


’  Casa  1:  Initialize  orbital  elements  for  time  of  epoch. 

’  Casa  2:  Update  orbital  elements  and  compute  in-plane  position  k  velocity. 
*  Simplest  of  NAVSPASUR  methods.  Accurate  to  10-16  n.mi.  within  20  days 
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of  epoch. 


*  «1m(1)  *  amajor 

'  elem(2)  a  ace 

’  elen(3)  «  tzaro 

’  elem(4)  »  aincl 

*  alam(S)  *  anode 

’  alam(G)  »  pari 

’  elem(7)  *  notion 

’  elem(8)  ■  manomaly 

’  elem(9)  *  decay 


a ami -major  azia 

eccentricity 

tine  of  perifocal  paaaage 

inclination  (radiana) 

longitude  of  aacending  node  (radiana) 

argument  of  perifocua  (radiana) 

naan  notion  (revolutiona/day) 

naan  anomaly  (radiana) 

let  decay  conatant  (radiana/minute*2) 


CONST  pi  s  3 . 1415926535897938 ,  rad  ■  pi  /  1808,  twopi  «  pi  ♦  pi 
CONST  c  «  .000811973858  *  0.76*j2 

CONST  twothird  «  28  /  38 

CONST  ae  s  6378.148  1  earth 'a  equatorial  radiua  (km.) 

CONST  GE  3  398600.58  ’  geo.  grav.  conat.  (km)‘3/(aec)*2 

k  «  608  e  SQR(GE  /  ae)  /  ae’ (eru)‘(3/2)/min 
CONST  we  3  1.0027379093507958  •  tuopi  /  14408’ aider eal  period  in  radiana/min 


SELECT  CASE  IX 
CASE  1 

elemt(l)  =  (k  /  elemi(7))  *  twothird 

a2#  3  elemt(l)  *  2 

ecc2t  =  1#  -  elem8(2)  *  2 

ci*  =  C0S(elem8(4)) 

c2i*  3  cif  «  ci# 

elem#(l)  =  elemt(l)  *  (18  ♦  e  *  (38  *  c2i8  -  18)  /  . 

(a28  *  ecc28  *  1.58))  *  twothird 
elem8(3)  3  time#  -  elem#(8)  /  elem#(7) 
aecc228  =  (elem#(l)  *  ecc28)  *  2 
peridot#  3  c  *  (58  *  c2i8  -  18)  /  aecc228 
nodedot#  =  -28  *  c  *  ci#  /  aecc228 
adot#  *  -28  3  twothird  *  elem#(l)  *  m28  /  elem#(7) 

FOR  iX  «  1  TO  9 

elemu#(iX)  *  elem#(iX) 

NEXT 
CASE  2 

ma.mO#  =  (elem#(9)  3  time#  *  elem#(7))  *  time# 

ma#  3  ma.mO#  ♦  elem#(8) 

elemu#(6)  =  elem#(6)  ♦  peridot#  *  ma.mO# 

elemu#(5)  3  elem#(5)  *  nodedot#  *  ma.mO#  -  we  *  time# 

elemu#(l)  =  elem#(l)  *  adot#  *  time# 

CALL  kepler(ma8,  elem#(2),  ea8) 
a08  3  SIN(ea8) 
c08  3  C0S(ea8) 
elemu#(8)  ■  ma# 

elemu#(3)  *  time#  -  ma#  /  elem#(7) 


’  in-plane  poaition  and  velocity 
aat.z#(l)  «  elemu#(l)  3  (cO#  -  elen#(2)) 
aat . z#(2)  «  eleau#(l)  3  SQR(ecc2#)  *  aO# 
aat.z#(3)  3  0# 

r#  ■  elemu#(l)  *  (1#  -  elem#(2)  3  c08) 

ed#  *  elemu#(l)  •  elem#(7)  /  r# 

aat.zdS(l)  *  -elemu#(l)  *  ed#  •  aO# 

aat.zd#(2)  •  eleau#(l)  *  ed#  *  SQR(ecc28)  *  cO# 

aat . zd#(3)  3  0# 

’  equatorial  rectangular  coordinatea 

CALL  euler(elemu(6)t  elemu(S),  elem(4)>  aat.zO,  sat.zdO) 

’  correct  velocity  for  earth ’a  rotation 
aat.zd#(l)  *  aat.zd#(l)  ♦  we  *  aat.z(2) 
aat.zd#(2)  ■  aat.zd#(2)  -  we  «  aat.z(l) 
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CASE  ELSE 
PRINT 

PRINT  "Error  in  pos. update.  First  argument  oust  ba  1  or  2  only." 
PRINT  "Argument  *  IX 
PRINT 
END  SELECT 
END  SUB 


SUB  risa.sat  (zv#(),  sta.RtO,  vat,  tO#,  aa#,  elenu#()) 

>  Compute  accantric  anomaly  for  tha  tima  of  risa  or  aat  of  a  satallita 
’  02-04-89.  Rav.  02-08-89  C  1700. 

’  NOTE:  To  nsa  this  routine,  tha  tiaa  of  culmination  MUST  ba  computed  first, 
'  Tha  antaring  aa#  argonaut  should  ba  aa.cnlt  -/♦  0.2  radians  (about  11 
’  degress)  for  risa/sat. 

’  See  the  introductory  comments  in  SUB  culminate  for  nomenclature  and 
'  notation. 

CONST  pi  =  3 . 141592653589793# ,  rad  »  pi  /  1808,  tsopi  »  pi  ♦  pi 

arnajt  =  elemu#(l) 
ecc#  -  elemu#(2) 
capt*  a  elemu*(3) 
incl*  *  «lemu#(4) 
node*  =  elemu#(5) 
perit  =  elemu#(6) 
mott  *  elemut(7) 

DIM  pv#(3) ,  pvdf(3),  qv# ( 3 ) ,  q»d#(3),  rvecd#(3) ,  rho#(3) 

cose#  «  C0S(peri#) 
sine#  a  SIN(parit) 
cosi#  a  COS(inclt) 
sini#  a  SIN(incl#) 

corr#  a  i 
DO 

cosea#  a  COS(eaS) 
ecosea#  =  ecc#  *  cosea# 
sinea#  a  SIN(ea#) 

esinea#  =  acc#  *  sinaa# 

gamma#  -  node#  -  set  *  ((aa#  -  esinea#)  /  mot#  ♦  capt#  -  t0#) 
gammad#  a  -ee#  *  (1#  -  ecosea#)  /  mot# 

cosg#  *  COS(gammaS) 
sing#  «  SIN(gamma#) 

rtacc#  «  amaj#  a  SQR(1#  -  acc#  a  acc#) 

xs#  ■  amaj#  a  (cosea#  -  acc#) 
xsd*  ■  -amaj#  «  sinaa# 
ys#  ■  rtacc#  *  sinaa# 
ysd#  «  rtacc #  *  cosaa# 

pv#(l)  «  cose#  a  cosg#  -  sins#  a  sing#  *  cosi# 
pv#(2)  *  cose#  a  sing#  ♦  sins*  a  cosg#  *  cosi# 
pv#(3)  *  sins#  a  sini# 

qv#(l)  '  -sins#  a  cosg#  -  coss#  •  sing#  *  cosi# 
qv#(2)  -  ins#  a  sing#  ♦  coss#  a  cosg#  •  cosi# 
qv#(3)  >  coss#  a  sini# 


pvd#(l)  »  -pv#(2)  *  gaamad# 
pvd#(2)  a  pv*(l)  *  gamaad# 
pvd#(3)  »  0# 

qvd#(l)  »  -qv#(2)  *  gaamad* 
qvd#(2)  *  qvf(l)  *  gaamad# 
qvd«(3)  *  0# 

FOR  iX  «  1  TO  3 

rho#(iX)  ■  xw»  *  pv#(iX)  ♦  )nl  *  qv#(iX)  -  sta.Rt(iX) 
rvecdS(iX)  ■  xwd#  *  pv#(iX)  ♦  xw#  *  pvd#(iX)  ♦  ywd#  *  qvS(iX)  _ 

♦  yw#  *  qvdt(iX) 

NEXT 

fa#  *  dot#(zv#(),  rho#()) 
fad#  »  dot#(zv#(),  rvacd#()) 

corr*  a  fa#  /  fad# 
aa#  a  «&#  -  corr# 

LOOP  UNTIL  ABS(corr#)  <  .00001# 

END  SUB 

> ****************************************************************************** 
SUB  sat. head  (r#,  sat.x#(),  sat.xdSO,  haad#)  ’  01-18-89.  Rav.  03-12-89. 

’satellite  heading 

CONST  pi  *  3 . 141592653589793# ,  tvopi  a  pi  ♦  pi 

ahead#  =  (sat.x#(l)  *  sat.xd#(2)  -  sat.x#(2)  *  aat.xd#Cl))  /  r# 
chead#  =  sat.xd#(3)  -  (sat.x#(3)  *  dot#(sat . x# ( ) ,  sat.xd#()))  /  r*  *  2 
haad#  a  atan2#( ahead#,  chead#) 

IF  haad#  <  0#  THEN  haad#  ■  head#  ♦  tvopi 

END  SUB 

>«**«•*********•***«•*«•*•*••*********••****••*•*•***•••*•**•*•*••*•«***••«**•* 
SUB  sat.  sphere,  coord  (aat.xO,  r,  t  j ,  tm,  sat.dc,  sat. In,  sat.ra) 

'compute  spherical  coordinates 

CONST  pi  *  3.141592653589793*,  rad  *  pi  /  180*,  tvopi  *  pi  ♦  pi 

sat.dc  a  asin(sat . x(3)  /  r)  ’  satellite  declination 

sat. In  a  atan2( sat . x(2) ,  sat.x(l))  '  satellite  longitude 

ru  a  iB«  «  gmst(tj,  ta)  •  rad  ’  Greenwich  mean  sidereal  time 

sat.ra  a  dmod(sat.ln  ♦  ru,  tvopi)  '  satellite  right  ascension 

END  SUB 


SUB  sta.p. coord  (It#,  In#,  ht#,  x#())  *  08-06-88.  Rev  02-04-89. 

’  convert  station,  lat,  long  fe  height  (ft.)  from  geographic 
'  to  geocentric  rectangular  -  x(.) 


CONST  nai .ft  a  .3048#  /  1852# 
CONST  km. ft  *  .3048#  /  1000# 
CONST  ae  «  6378.14# 

CONST  fl  >  i#  /  298.257# 

CONST  e.sqr  ■  (2#  -  fl)  a  fl 


*  a. mi. /foot 
’  km. /foot 

’  earth’s  equatorial  radius  (km.) 
’  earth’s  flattening  factor 

*  earth’s  allipticity  squared 


ht.eru#  ■  ht#  •  km. ft  /  ae 

sit#  «  SIN(lt#) 

clt#  «  C0S(lt#) 

sin#  >  SIN(ln#) 

cln#  ■  C0S(ln#) 

n#  ■  1#  /  SQR(1#  -  e.sqr  •  sit  *  sit) 

gl#  ■  n*  ♦  ht.eru# 


40 


g2 #  *  nt  *  (It  -  • . »qr)  ♦  ht . eru# 


'radius  vector  lor  flattened  earth 
*i(l,  1)  »  gl*  *  cltf  *  cln# 

x#(2,  1)  *  gl*  *  clt*  *  a In# 

x*(3,  1)  *  g2t  •  alt* 

’ Spherical  aarth  unit  vactora: 

’l)  local  parpandicular  vactor 
x#(l,  2)  «  clt#  *  cln# 

x#(2,  2)  >  clt#  •  aln# 

x#(3 ,  2)  *  alt# 

*2)  aaat  vactor 

x#(l,  3)  =  -aln# 

x#(2,  3)  a  cln# 

x#(3,  3)  »  0# 

’3)  north  vactor 
x#(l,  4)  =  -*lt#  *  cln# 

x#(2,  4)  *  -alt#  *  aln# 

x#(3,  4)  =  clt# 

END  SUB 

>*•••«*>«•«**«*«•«••****•*****•**••*••••***««•*•««*••*•••**••**«******••*••*•••• 
SUB  sta.sat  (ata.x#(),  aat.x#(),  rngt,  az# ,  el#,  langla#) 

■  08-06-88.  Rav.  08-31-88. 

’compute  atation-eatellita  ralationa.  Range,  azimuth,  elevation  and  satellite 
’  "look-angle"  (angle  batveen  aatallite  aub  point  and  atation) . 

CONST  pi  a  3.141592653589793#,  rad  *  pi  /  180#,  twopi  *  pi  ♦  pi 

’rho()  a  sat  coords  ainua  ata  coorda.  rhodota  are  the  rho  components  projected 
’  in  the  station  local  north-aaat-radial  system, 

DIN  rho#(3),  rhodots#(3) 

mg#  =  0* 

FOR  iX  =  1  TO  3 

r#  a  sat.x#(i%)  -  sta.x#(i%,  1) 

rho#(iX)  =  r* 

rng#  *  rng*  +  r»  *  r# 

NEXT 

rng#  »  SQR(rng#) 

FOR  iX  «  1  TO  3 

rho#(iX)  a  rho#(iX)  /  rng# 

NEXT 

FOR  iX  »  1  TO  3 

a#  a  0# 

FOR  jX  »  1  TO  3 

a#  «  a#  ♦  rho#(jX)  *  sta.x#(jX,  it  ♦  1) 

NEXT 

rhodota#(iX)  *  a# 

NEXT 

az#  a  atan2#(rhodots#(2) ,  rhodota#(3))  /  rad 
IF  az#  <  0*  THEN  az*  ■  az#  ♦  360# 

a#  ■  rhodots#(l) 

al*  a  ATN(a#  /  SQR(1#  -  a#  *  a#))  /  rad 

a#  a  0# 
b#  a  0# 

FOR  iX  ■  1  TO  3 
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c t  *  sat.xt(iX) 

at  »  at  ♦  rhot(iX)  *  c t 

bt  »  bt  ♦  ct  *  ct 

HEXT 

langlat  ■  acoat(at  /  SQR(bt))  /  rad 
END  SUB 


SUB  tima.npdata  (kX,  dalta.tt,  tat,  jdt,  jdt(),  tjt,  tiaa.of.day,  at,  dt,  _ 
yk,  ht,  ao$,  day!) 

CONST  twantyfour  *  24 t,  aixty  *  60t,  jul.csnt  ■  36625# 

t#  a  dalta.tt  'dalta. tt  in  hour*  if  kX  ■  0,  tt  in  hours. 

IF  kX  >  0  THEN  tt  *  tt  /  sixty  'dalta.tt  in  ainntas  if  kX  >  0,  tt  in  hours, 
tat  *  tat  ♦  tt 

jdt(2)  «  jdt(2)  ♦  tt  /  twantyfour 

jdt  *  jdt(l)  ♦  jdt(2) 

tjt  -  tjt  ♦  tt  /  twantyfour  /  jul.cant 

tima.of.day  a  tima.of.day  ♦  tt 

IF  tima.of.day  >=  twantyfour  THEN 

tima.of.day  a  tima.of.day  -  twantyfour 

jdt(l)  =  jdt(l)  ♦  It:  jdt(2)  »  jdt(2)  -  It 

CALL  jd.to.data( jdt() ,  mk,  dk,  yk ,  ht,  mot,  dayt) 

END  IF 

IF  tima.of.day  <  Ot  THEN 

tima.of.day  -  tima.of.day  ♦  twantyfour 

jdt(l)  a  jdt(l)  -  It:  jdt(2)  =  jdt(2)  ♦  it 

CALL  jd.to.dataCjdtO ,  mk,  dk,  yk,  ht,  mot,  dayt) 

END  IF 

END  SUB 

'»•«««**»«****•***********•**«***««******•**•••••**•«*•*•**•**%***•*******••**« 
SUB  yrmoday  (vt,  yrk,  mot,  dyk)  '  03-21-88.  Rav  03-21-88. 

’  Convart  yyyy.mmdd  to  yrk=yyyy,  mok=mm,  and  dyk*dd 

IF  INSTRCvt.  <>  5  OR  LEN(vt)  <>  9  THEN 

PRINT 

PRINT  "Error  in  yrmoday.  vt  must  ba  of  tha  form  yyyy.mmdd" 

PRINT  "vt  =  ";  vt 
PRINT 

ELSE 

yrk  =  VAL(LEFTt(vt,  4)) 
mok  a  VAL(HIDt(vt,  6,  2)) 
dyk  a  VAL(RIGHTt(vt,  2)) 

END  IF 
END  SUB 


APPENDIX  E.  THE  HIGH  ACCURACY  KINEMATICS  ROUTINE. 

This  appendix  contains  a  listing  of  the  alternate  pos. update  routine.  It  is  a  BASIC  ver¬ 
sion  of  the  highly  accurate  PPT2  subroutine  provided  by  NAVSPASUR  [Ref.  1]  with  their 
SHOWALL  program.  It  is  believed  that  it  has  been  accurately  translated,  however  any  errors 
are  solely  those  of  the  author.  Known  discrepacies  are  the  difference  in  the  WGS-72  con¬ 
stants  used  in  SHOWALL  and  the  IAU  1976  constants  used  here  [see  section  III.  SOURCES, 
in  this  report  for  details]. 


SUB  pos. update  (indX,  elem9(),  tin**,  sat.x9(),  sat.xd9(),  elemu9(),  _ 
rt,  •«.#)  STATIC  Rev.  09-26-89  «  0900 

’  Casa  1:  Initialize  orbital  elements  for  time  of  epoch. 

'  Casa  2:  Update  orbital  alamanta  and  compute  in-plane  position  k  velocity. 

’  Comments  from  original  FORTRAN  version: 

*  implicit  real*8(a-h ,o-z) 

1 ***  model  ref erence-Astronomical  Journal  64, No  1274, Nov. , ’59 
•••*  solution  of  the  problem  of  artificial  satellite  theory 
’***  with  out  drag... Dirk  Brouwer. 

'  840329  confirmed  use  of  w*u  x  v  vice  »*n  s  vel 

'  840329  inserted  cdh.sdh,  cdt.sdt  re  osc  vice  mean  els  at  t 

’  for  use  only  with  dcext2,dcl,dcia 

’  frame  for  u,v,w  is  the  rotated  inertial  w/o  coriolis 
'  830127  new  partials  re  (a,x,y ,z,p,q) 

’  830127  improved  (polynomial)  model  for  singularity  in  Brouwer 

DIN  a(25) 

DIN  f (25) ,  osc(lO),  kf (10) ,  cf(10),  bs(3,  4),  u(3),  v(3),  w(3),  vel(3) 

DIN  ar(6 ,  8),  br(3,  6),  cr(3,  6),  dv(3) 


’**»  The  commented  beta  value  is  the  WGS-72  value .however  using  this 
>«**  value,  one  obtains  a  value  for  a(9)  different  then  what  we  had 
'*•*  been  using,  thus  causing  a  problem  with  tim» .therefore  to  retain 
’***  previous  value  of  a(9),one  solves  for  beta  using  the  old  value 
'***  of  a(9) . 

>*«*  betal*398600 . SdO 

’  betal  *  398597.625795889,  erkm  »  6378.1359,  flat  *  298.269 
*  k  »  .07436691619 


'  IAU  1976: 

CONST  twopi  -  6.2831853071795869 
CONST  we  *  1.0027379093507959  *  twopi  / 
CONST  ae  *  6378.149,  erkm  »  ae 
CONST  betal  *  398600.59 

k  *  609  *  SQR(betal  /  ae)  /  ae 
CONST  flat  «  19  /  298.2579 


14409’sidereal  period  in  radians/min 
’  earth's  equatorial  radios  (km.) 

'  geo.  grav.  const.  (km)*3/(sec) *2 

’(eru)*(3/2)/min 

'  earth's  flattening  factor 


’  k-terms  [Note:  These  may  or  may  not  agree  with  IAU  1976] 

CONST  c20  >  -.00048416059,  c30  ■  .000000959589 
CONST  c40  ■  .000000551999,  c50  «  .0000000658759 


CONST  twotrd  >  .6666666666679,  fortrd  «  1.3333333333339 
CONST  tentrd  ■  3.3333333333339,  ar  ■  48  •  09 


'IF  indX  <>  1  THEN  GOTO  secondcase 
SELECT  CASE  indX 
CASE  1 

kzX  «  1 

a(l)  *  09'  [Note:  ■  85-01-01  in  SHOWALL.  Set  to  zero  here.] 
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a(3)  *  -.8#  *  c20  *  SQR(8#) 
a(4)  «  c30  *  SQR(7#) 
a(5)  »  .375  *  c40  •  SQR(9#) 
a(6)  «  c50  *  SQR(llt) 

’  two  pi 
a(7)  «  twopi 
a(l6)  «  1#  /  a(7) 

’  hergs/day  ,  secs/herg,  mins/herg 
a(9)  »  arks  *  SQR(erkm  /  betal) 
a(8)  *  86400*  /  a(9) 
a(l7)  *  1440#  /  a(8) 

’  we(rad/day),  we-  2  pi,  we(rad/herg) 
a(ll)  *  6.30038809878 
a(10)  =  a(ll)  -  a(7) 
a(l2)  *  a(ll)  /  a(8) 

’  aarth  flattening 

a(13)  «  (2#  -  1#  /  flat)  /  flat 

’  sm/er,  ka/ar,  nm/er 

*(20)  =  arkm  /  1.609344# 

a(21)  *  arka 

a(22)  =  arka  /  1.852# 

’  dag/rad 
a(23)  *  360#  /  twopi 

range  rate/er/herg  to  cyclaa/aacond  -  conversion 
a(24)  =  a(21)  *  216980000#  /  (a(9)  *  299792.5#) 

’  fanca  plana  displac  from  aarth  cantar 
a(25)  =  .0031# 

’  Patch  to  interface  PPT2  and  SATSTA  variable*: 

f(l)  =  alem(8)  'mean  anomaly 

f (2)  =  alam(7)  a  a(l7)  'mean  motion 

f(3)  =  alam(9)  *  a(17)  *  2’firat  decay  conatant 

f(4)  =  0#  ’second  decay  conatant 

f(5)  =  elam(2)  'eccentricity 

f(6)  -  alam(6)  'arg  of  perigee 

f(7)  *  elem(5)  ’node 

f(8)  =  C0S(alem(4))  ’cosine  inclination 

f(9)  ~  0#  ’  epoch  in  PPT2,  time  from  epoch  hare  (?  but  check  ?) 

’  Used  only  once  to  "recover"  the  remaining  elements 
hi  =  0# 

'***  f(8)  is  cosine  of  inclination 

t9  =  f (8)  *  f (8) 
tl  =  1#  -  5#  *  t9 

’  t2=(1.0d0-axp(-100.0d0*tl*2))/tl 

revised  handling  of  divisor  par  navspasur  o2t  memo  nov  '83 
beta  *  100#  /  2#  *  11 
p3  «  beta  *  tl  *  2 
pi  «  EXP(-p3) 
p2  *  1#  ♦  pi 
pier  «  pi  *  2 
p4  *  1# 

p3ex  *  -.8#  *  p3 

FOR  nX  *  2  TO  13 

IF  nX  <■  11  THEM  p2  «  p2  *  (1#  a  pl*x) 

plax  •  plax  *  plax 

p4  ■  p4  ♦  p3ax 

p3ax  ■  -p3  *  p3ax  /  (nX  ♦  1) 

REIT 

t2  *  p2  *  p4  *  beta  *  tl 

cio2  ■  SQR( . 6#  ♦  .6#  *  f(8)) 

sio2  «  SQRC.5S  -  .8#  a  f(8)) 

tl2  »  3#  *  t9  -  1# 

tl3  •  t9  a  t9 

tl4  »  t2  a  tl3 

tl5  »  1#  -  t9 
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'***  f(5)  is  eccentricity 

•sqO  *  1(5)  *  1(5) 

et*20  *  it  -  ssqO 

etaO  *  SQR(eta20) 

t6  »  staO  *  sta20 

t7  «  ata20  *  ata20 

t8  ■  ataO  ♦  it  /  (it  ♦  ataO) 

1(13)  *  1(2)  *  (-twotrd) 

•ls«(l)  -  1(13) 

'***  1(2)  is  basically  nean  notion  but  also  contains  tha 

’***  sscular  term  in  Bean  anoaaly.  therefore  1(13)  is  (alaost)  the 

’***  seai  major  axis. 

’***  This  loop  is  to  reaove  the  contribution  ol  the  additional 
’***  term  in  aean  aotion  lor  iaproving  values  ol  seai-Bajor  axis 
'***  and  seculars. 

’***  a(3)  to  a(6)  are  ’k’  terms  wgs  72 

FOR  iX  ■  1  TO  6 

lsl3  «  1(13)  *  1(13)  *  t7 

Is  1 1  a  a(3)  /  lsl3 

Is  13  *  a(5)  /  (lsl3  *  isl3) 

tt  *  It  ♦  1.5t  *  Isll  *  etaO  *  tl2 

uu  =  -15t  ♦  16f  •  etaO  ♦  25f  •  eta20 

uu  =  uu  ♦  (30t  -  96t  *  etaO  -  90t  *  eta20)  *  t9 

uu  ■  uu  ♦  (105*  ♦  144*  *  etaO  ♦  25*  *  eta20)  *  tl3 

tt  =  tt  ♦  .09375*  *  Isll  *  2  *  etaO  *  uu 

tt  =  tt  ♦  .9375*  *  lsl3  *  etaO  *  esqO  *  (3*  -  30*  *  t9  ♦  35*  *  tl3) 

1(13)  =  (tt  /  1(2))  *  twotrd 

NEXT 

’***  1(11)  sill  be  rate  ol  change  ol  arguaent  ol  perigee 

uu  -  -35*  +  24*  *  etaO  ♦  25*  *  eta20 

uu  -  uu  ♦  (90*  -  192*  *  etaO  -  126*  *  eta20)  *  t9 

uu  =  uu  ♦  (385*  ♦  360*  *  etaO  +  45*  *  eta20)  *  tl3 

tt  «  -1.5*  *  Isll  *  tl  ♦  .09375*  *  Isll  *  2  *  uu 

uu  =  21*  -  9*  *  eta20  ♦  (-270*  ♦  126*  *  eta20)  *  t9 

uu  *  uu  ♦  (385*  -  189*  *  eta20)  *  tl3 

f(ll)  »  tt  ♦  .3125*  *  1  sl3  *  uu 

’***  1(12)  will  be  rate  ol  change  ol  right  ascention 

uu  =  -5*  +  12*  *  etaO  ♦  9*  *  eta20 

uu  =  uu  -  (35*  ♦  36*  *  etaO  ♦  5*  *  eta20)  *  t9 

tt  =  -3*  *  Isll  *  1(8)  ♦  .375*  *  Isll  "  2  *  1(8)  *  uu 

1(12)  ■  tt  ♦  1.25*  a  1 sl3  *  1(8)  *  (6*  -  3*  *  eta20)  *  (3*  -  7*  *  t9) 

’***  l(15)=sine  inclination, l(l4)is  a  dot, 1(16)  is  eccen.dot 
1(15)  =  SQR(t 15) 

1(14)  *  -lortrd  *  1(3)  a  1(13)  /  1(2) 

1(16)  »  1(14)  *  1(5)  *  eta20  /  1(13) 
ql  »  1(2)  *  1(13)  *  1.5* 

1(11)  *  1(11)  /  ql 

1(12)  «  1(12)  /  ql 

lsl2  «  a(4)  /  (a(3)  *  1(13)  *  eta20) 

lsl3  •  lsl3  /  Isll  •  tentrd 

1 sl4  ■  a(6)  /  (a(3)  *  1(13)  *  3  *  eta20  *  t7) 

ql  «  .126*  *  (Isll  *  (1*  -  11*  •  t9  -  40*  *  tl4)  -  lsl3  . 

*  (1*  -  3*  *  t9  -  8*  •  tl4)) 

p5  «  1*  ♦  t2  •  (8*  *  t9  ♦  20*  •  t 14) 

p2  »  1*  ♦  2*  *  pS 

q2  »  .126*  •  esqO  •  1(8)  *  (Isll  *  (1*  ♦  10*  •  p5)  -  1 sl3  *  p2) 

p2  >  .46875*  *  p2  •  1(5)  •  1(8)  «  1(16)  *  (4*  ♦  3*  *  esqO)  •  lsl4 

p3  >  lsl4  *  (1*  -  9*  *  t9  -  24*  *  tl4) 
q6  a  .25*  e  (lsl2  ♦  .3125*  a  (4*  ♦  3*  *  esqO)  •  p3) 
p3  ■  .15625*  *  1(5)  *  1(15)  a  p3 

p4  a  .030381944*  a  1(5)  a  1,14  a  (1*  -  5*  a  t9  -  16*  a  tl4) 

p6  a  .060763889*  a  1(5)  a  esqO  a  1(8)  a  1(15)  a  lsl4  a  (1*  ♦  4*  a  pS) 

vlel  ■  1(6)  *  eta20  •  ql 

vlhli  ■  -1(15)  a  q2 

tt  »  (t6  -  1*)  *  ql  -  q2  ♦  26*  *  esqO  a  tl4  •  t9  a  t2  •  (Isll  -  .2*  •  lsl3) 
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vlsl  «  tt  -  .0625#  *  asqO  *  (f*ll  *  (1#  -  33#  *  t9  -  200#  *  tl4) 

-  Isl3  *  (1#  -  9#  *  t9  -  40#  *  tl4)) 
vl*2  *  ata20  *  1(15)  *  q5 
vlh2i  »  1(5)  *  1(8)  *  q5  ♦  1(15)  *  p2 
vl«2  «  1(5)  *  1(16)  *  (t8  ♦  1(8)  /  (1#  ♦  1(8)))  *  qS  _ 

♦  (11#  ♦  3#  *  asqO  -  3#  *  t6)  *  p3  ♦  (l#  -  1(8))  *  p2 

vla3  =  -3#  *  1(5)  *  ata20  *  1(16)  *  p4 
vlh3i  «  -asqO  *  1(8)  *  p4  -  1(16)  *  p6 

vl«3  *  1(15)  *  (3#  *  t6  -  3#  -  2#  *  asqO  -  asqO  *  1(8)  /  (1#  ♦  1(8)))  *  p4 

-  (1#  -  1(8))  *  p6  ^ 

vll2  *  vla2  ♦  3#  *  1(5)  *  ata20  *  p3 
'  pracomps  lor  tha  partial* 

’***  *(6)  ia  argnaant  ol  parigaa,l(7)  ia  right  aacanaion 

ainhO  «  S1H(1(7)) 

coahO  *  C0S(1 (7)) 

aintO  »  SI*(1(6)  ♦  1(7)) 

coatO  «  C0S(1(6)  ♦  1(7)) 

’  kaplar’a  law 

dada  «  -tsotrd  *  1(13)  /  1(2) 

1  variation*  thru  th*  aaculara 

aacul  =  a(3)  /  (1(13)  *  wta20)  *  2 

dgddm  *  3#  *  aacul  *  tl  /  1(13)  *  dada 

dgdde  =  -6#  *  aacul  *  tl  *  1(5)  /  *ta20 

dgddi  a  -is#  *  a«cul  *  1(8)  *  1(15) 

dhdda  a  6#  *  aacul  *  1(8)  /  1(13)  *  dada 

dhdda  a  -12#  *  aacul  *  1(8)  •  1(5)  /  ata20 

dhddi  a  3#  *  aacul  *  1(15) 

dtddz  a  (dgdda  +  dhdda)  •  coatO 

dtddy  a  (dgdda  ♦  dhdda)  *  aintO 

dtddp  a  (dgddi  ♦  dhddi)  *  2#  *  coahO  /  cio2 

dtddq  a  (dgddi  ♦  dhddi)  *  2#  *  ainhO  /  cio2 

dhddz  a  dhdda  «  coatO 

dhddy  a  dhdda  *  aintO 

dhddp  a  dhddi  *  2#  *  coahO  /  cio2 

dhddq  a  dhddi  *  2#  *  ainhO  /  cio2 

’  variation*  thru  m2 

daddml  a  lortrd  *  1(3)  *  (1(13)  /  1(2)  -  dadm)  /  1(2) 

daddm2  a  -lortrd  *  1(13)  /  1(2) 

dwd<i*  a  -lortrd  »  (1#  -  3#  »  aaqO)  *  1(3)  /  1(2) 

daddz  a  dadda  •  costO 

daddy  a  dadda  *  aintO 

daddml  =  lortrd  *  1(5)  *  ata20  *  1(3)  /  1(2)  *  2 
daddm2  a  -lortrd  *  1(5)  *  ata20  /  1(2) 

FOR  iX  =  1  TO  9 

alamu#(iX)  a  *l*m#(iX) 

HEXT 

'aacondcaaa: 

CASE  ELSE 

tm  »  tima#  /  u(17)  'Convert  minuta*  to  Bargs 
hi  «  ta  -  1(9) 

'***  computa  position  as  a  Inaction  ol  tima  (in  call  lina) 
osc(9)  a  hi  •  (1(2)  ♦  hi  *  (1(3)  ♦  hi  a  1(4))) 
ose(3)  a  dmod(l(l)  ♦  oac(9),  twopi) 
tt  a  1(5)  ♦  1(16)  a  hi 
IF  tt  <  0#  THE*  tt  a  0# 

IF  tt  >  .99999#  THEN  tt  «  .99999# 

oac(4)  »  tt 

asq  »  osc(4)  *  2 

ata2  *  1#  -  asq 

ata  «  SQR(ata2) 

CALL  kaplar(osc(3) ,  o*c(4),  e3) 
osc(l)  a  C0S(c3) 
cl  *  1#  -  osc(4)  •  osc(l) 
osc(l)  «  (osc(l)  -  osc(4) )  /  cl 
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oac(2)  =  «tt  *  SIN(c3)  /  cl 

ta4  *  1#  ♦  oac(4)  *  osc(l) 
cf 3  =  0# 

IF  kz%  <>  0  THEN  cf3  a  a(l2) 
tt  *  f  (13)  ♦  f (14)  *  hi 
IF  tt  <  1#  THEN  tt  *  1# 
oae(8)  =  tt 
oac(7)  a  f(8) 
fai  *  f  (15) 

obc(6)  «  dmod(f(7)  ♦  f(12)  *  oac(9),  t*opi) 
ogc(S)  *  dmod(f (6)  ♦  f(ll)  *  oac(9),  twopi) 

r  =  oac(8)  *  ata2  /  ta4 

obc(6)  =  oac(6)  -  cf3  *  (tm  -  »(1)) 


*1  = 

C0S(oac(5)) 

*2  = 

SIN(oac(5)) 

*7  a 

C0S(oac(6)) 

*8  a 

SIN(oac(6)) 

agda 

=  0# 

agda 

a  0# 

agdi 

=  0# 

agdl 

=  0# 

*gdg 

=  0* 

agdh 

=  0# 

*3  a 

*1  *  *1  -  *2  * 

*2 

*4  a 

2*  a  *1  *  *2 

*5  a 

*3  *  *1  -  *4  * 

*2 

*6  a 

*4  *  *1  *  *3  * 

*2 

*9  a 

oac(l)  *  08c( 1 ) 

-  osc(2)  a  oac(2) 

*10  a 

2#  *  os  c ( 1 )  * 

oac(2) 

*11  = 

*3  *  osc(l)  - 

*4  *  obc(2) 

*12  a 

*4  *  osc(l)  * 

*3  *  oac(2) 

*13  = 

*3  *  *9  -  *4  * 

*10 

*14  *  »4  *  »9  ♦  »3  *  *10 

*15  s  *13  *  oac(l)  -  *14  *  oac(2) 

*16  *  *14  *  oac(l)  ♦  *13  *  oac(2) 

*17  =  atan2(osc(2),  oac(l))  ♦  oec(4)  *  oac(2)  -  dmod(oBc(3) ,  t*opi) 

*18  =  C0S(oac(3) ) 

*19  =  SIN(oac(3)) 

*20  a  oac(l)  *  (3#  ♦  oac(4)  *  oac(l)  *  (3#  ♦  oac(4)  *  oac(l))) 

*21  =  3#  *  *14  +  3*  *  oac(4)  *  *12  ♦  oac(4)  *  *16 
*22  =  ta4  *  (1#  ♦  t s4)  /  ata20 

oac(8)  =  oac(8)  *  (1#  +  fall  /  *ta20  *  (tl2  *  (ta4  ‘  3  -  t6) 

♦  3*  *  tlS  «  ta4  ‘  3  *  *13))  ♦  agda 
da  =  vl*l  *  *3  ♦  vl*2  *  *2  ♦  vl*3  *  *6 

di  a  aio2  ♦  .5#  *  cio2  *  (.5#  *  fall  *  f(8)  *  f(15)  *  (3#  *  *13  +  3*  *  f(5) 

*  *11  ♦  f(S)  *  *1S)  -  d*  *  f(S)  *  f (8)  /  f (15)  /  ata20)  ♦  .68  *  cio2  *  agdi 

d*  a  oac (4)  ♦  .5#  *  fall  *  (tl2  *  (*20  ♦  f(5)  *  t8)  ♦  3#  a  tlS  *  (*20  ♦  f(5)) 

*  *13  -  *ta20  a  tlS  *  (3#  •  *11  ♦  *16))  ♦  da  ♦  agda 

dh  a  .5*  /  cio2  *  (-.5#  *  fall  a  f(8)  a  f(lB)  a  (6*  a  *17  -  *21)  ♦  vlhli  a  *4 

♦  vlh2i  *  *1  ♦  vlh3i  a  *5)  t  . 5f  a  agdh  /  cio2 
dl  a  -.25#  a  t6  a  fall  a  (2#  a  t12  a  (,22  ♦  1#)  a  0ac(2)  ♦  3#  a  tl5  •  ((1# 

-  *22)  a  *12  ♦  (.333333333#  ♦  *22)  a  *16)) 
oa  *  oac(3)  ♦  oac(5)  ♦  oac(6)  -  dl  a  f(8)  a  (t8  -  1#)  /  t6  -  .25#  a  fall  a  (6# 

a  *17  a  (tl  ♦  2#  a  f(8))  -  *21  a  (tl  ♦  2#  ♦  2#  *  f(8)))  . 

♦  vial  a  si  ♦  *i,2  a  *1  ♦  *1,3  a  ,5  ♦  agdg 

dl  *  dl  ♦  ataO  a  (*l,i  a  *4  -  *H2  a  *1  -  *1,3  a  *5)  ♦  agdl 
•aq  ■  da  •  da  *  dl  *  dl 
oac(4)  ■  SQR(aaq) 

oac(3)  »  atan2(da  »  *19  ♦  dl  a  *18,  da  a  *18  -  dl  a  *19) 

oac(7)  a  if  -  2#  a  (di  a  di  ♦  dh  *  dh) 

fai  «  S2R(1#  -  oac(7)  a  oac(7)) 

oac(6)  «  atan2(di  *  *8  ♦  dh  *  *7,  di  a  *7  -  dh  •  *8) 

oac(S)  ■  oa  -  oac(3)  -  oac(6) 
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eta2  *  It  -  esq 

eta  a  SQR(eta2) 

CALL  kepler(oac(3) ,  oac(4),  c3) 
osc(l)  «  C0S(c3) 
cl  =  If  -  osc(4)  *  osc(l) 
osc(l)  a  (oac(l)  -  osc(4))  /  Cl 
osc(2)  «  at a  a  SIN(c3)  /  cl 

ts4  =  It  ♦  oac(4)  *  oac(l) 
r  =  oac(8)  «  ata2  /  ta4 
•7  a  C0S(oec(6)) 

*8  a  SIN(oac(6)) 
el  a  C0S(oac(5)) 
e2  a  SIN(osc(5)) 
uba  *  »2  *  oac(l)  ♦  el  *  oac(2) 
ubc  *  el  *  osc(l)  -  w2  *  oac(2) 
n(l)  a  »7  *  nbc  -  e8  *  aba  *  oac(7) 
u(2)  «  *8  a  ubc  ♦  e7  *  aba  *  oec(7) 
u(3)  *  aba  *  fai 

v(l)  a  -*7  a  uba  -  w8  *  ubc  *  oac(7) 

v(2)  =  -w8  *  uba  ♦  *7  *  ubc  *  oac(7) 

v(3)  =  ubc  *  fai 

»(1)  =  v8  *  fai 

w(2)  =  ~w7  a  fai 

b(3)  *  oac(7) 

ta2  a  #ta  *  SqR(oac(8)) 

FOR  iX  *  l  TO  3 

vel(iX)  =  (osc(4)  *  oac(2)  »  u(iX)  ♦  ta4  *  v(iX))  /  ta2 

NEXT 

vel(l)  »  val(l)  ♦  cf3  *  r  *  u(2) 
v®l(2)  a  v*l(2)  -  cf 3  a  r  *  u(l) 

*  Patch  for  PPT2  to  SATSTA  output 
«l«mu(l)  a  oac(8) ’semimajor  aria 
elemu(2)  a  oac(4) ' eccentricity 

’ elemu(3)a  ’time  of  last  perigee  passage 

elemu(4)  =  acos(oec(7)) ’inclination 

elemu(5)  a  oac(6)’node 

elemu(6)  a  oac(5)’arg  of  perigee 

elemu(8)  a  oac(3)’mean  anomaly 

mat  a  el emu (8) 

CALL  kepler(maf,  elemi(2),  eat) 
sOt  a  SIN(eat) 
cOt  a  COS(eaf) 

’elemut(8)  *  mat 

elemut(3)  *  timet  -  mat  /  elemf(7) 

’  in-plane  position  and  velocity 
ecc2t  a  it  -  elemut(2)  *  2 
sat.zf(l)  a  elemut(l)  *  (cOt  -  elemt(2)) 
aat.xt(2)  >  elemut(l)  a  SQR(ecc2t)  •  aOt 
sat .zf (3)  •  Of 

rt  a  elemut(l)  a  (l#  -  elemf(2)  a  cOt) 

edt  «  elemnt(l)  *  elemt(7)  /  rf 

aat.zdf(l)  a  -elemut(l)  a  edt  a  aOt 

sat . zdt(2)  •  elemut(l)  a  edt  a  SQR(ecc2f)  a  cOt 

sat.xdt(3)  ■  Of 

'  equatorial  rectangular  coordinates 

CALL  euler(elemu(6) ,  elemu(S),  elem(4),  aat.z(),  aat.zdO) 

END  SELECT 
END  SUB 
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